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ABSTRACT 


Factorization  is  an  approach  to  linear  programming  (LP)  in  which  the  algebraic 
elements  of  the  LP  tableau  are  organized  in  such  a  way  that  a  large  portion  of  the 
tableau  may  be  represented  implicitly  and  generated  from  the  remaining  explicit 
pari.  In  dynamic  row  factorization,  the  row  structure  of  the  LP  model  instance 
influences  the  algebraic  structure  of  the  tableau,  and  the  dimension  of  the  algebraic 
elements  may  change  as  the  solution  progresses. 

We  present  three  algorithms  motivated  by  this  approach,  each  resulting  from 
a  different  LP  model  row  structure;  generalized  upper  bound  (GUB)  rows,  pure  net¬ 
work  rows  and  generalized  network  rows.  We  describe  implementations  of  all  three 
algorithms,  specifying  data  structures  for  tableau  and  basis  inverse  representations 
and  detailing  procedures  for  manipulation  and  update  of  these  representations. 

Computational  results  are  presented  for  a  number  of  real-world  models  taken 
from  a  variety  of  applications  and  industries.  From  each  model,  one  or  more  partic¬ 
ular  instances  are  solved  by  each  of  otir  three  implementations  and  by  a  commercial- 
quality  mathematical  programming  system.  The  characteristics  of  the  four  solvers 
are  compared  and  contrasted.  Previous  research  on  related  algorithms  by  others 
suggests  that  tliese  algorithms  are  properly  viewed  as  specialized  approaches,  useful 
only  on  narrow  cla.^se.'-  of  problems.  Our  computational  results  strongly  refute  this 
view,  and  instead  suggest  that  each  algorithm  is  superior  to  the  general  simplex 
approach  on  a  w  idc  range  of  problem  classes  and  structures. 
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I.  INTRODUCTION 


A.  INTRODUCTION 

A  recurring  theme  in  the  development  of  simplex-based  algorithms  for  hnear 
programming  has  been  the  identification  and  exploitation  of  special  problem  struc¬ 
ture.  Ideas  as  apparently  disparate  as  the  simplex  method  for  bounded  variables, 
primal  and  dual  decomposition  methods,  pure  and  generalized  network  primal  sim¬ 
plex  algorithms  and  primal  partitioning  schemes  may  be  unified  to  a  degree  by 
interpreting  their  development  in  this  context. 

The  factorization  approach  due  to  Graves  and  McBride  [1976]  provides  a  unify¬ 
ing  framework  from  which  we  may  reinterpret  many  existing  algorithms  and,  through 
the  application  of  the  common  principles  embodied  in  the  approach,  develop  new 
algorithms.  This  approach  has  as  its  central  thesis  the  idea  of  recognizing,  isolat¬ 
ing  and  exploiting  sjX'cial  structures  which  may  occur  in  a  certain  type  of  linear 
programming  tableau.  .An  exainpb-  of  such  special  structure  is  generalized  network 
rows.  Cieneralizecl  network  rows  ar<.’  naturally  specified  as  a  row  structure  in  which 
each  column  lias  at  most  two  nonzero  elements  within  the  rows. 

This  paper  adojits  the  design  principles  and  algorithmic  structure  suggested 
by  Graves  and  McBride  [197(i  to  develop  algorithms  that,  while  general  in  nature  in 
the  sense  that  each  may  be  used  to  solve  any  linear  programming  problem  instance, 
are  strongly  tailored  td  e.xploil  a  particular  row  structure.  W'e  call  this  approach 
"d}uiamic  row  fact oi  i/ai  iiu.  :  "row  factorization  '  because  we  exploit  the  structure 
in  the  basic  tableau  v.  hn  I.  induced  by  the  row  structure  of  the  LP  model  instance, 
and  "dynamic"  because  t  he  dimension  of  the  structure  may  vary  (or  even  fail  to  be 


present)  as  the  solution  progresses.  In  our  setting,  we  require  the  row  structure  of 
the  model  instance  to  be  specified  prior  to  solving,  and  to  remaun  fixed  throughout 
the  solution  process.  An  extension  of  this  approach  is  to  allow  the  row  structure  of 
the  model  instance  to  vary  as  the  problem  is  solved.  While  we  do  not  develop  an 
implementation  of  the  second  approach,  we  show  that  it  is  a  conceptuadly  simple 
extension  of  our  work. 

Each  algorithm  is  developed  by  factoring  the  constraints  of  the  LP  model 
into  two  classes:  those  that  have  special  structure  (factored)  and  those  that  do 
not  (explicit).  This  factoring  of  constraints  induces  a  factored  structure  in  the  LP 
tableaus  which  may  be  exploited  computationally.  We  treat  each  of  three  structures 
of  factored  constraints  in  this  work:  generalized  upper  bounds  (GUB),  pure  networks 
and  generalized  networks. 

We  implement  each  of  the  three  algorithms  by  integrating  it  within  the  struc¬ 
ture  of  a  state-of-the-art  mathematical  programming  system:  the  X-System  of 
Brown  and  Graves  [1975j.  We  do  so  both  to  demonstrate  the  feasibility  of  imple¬ 
menting  such  methods  and  to  determine  whether  or  not  such  met  hods  have  practical 
value  for  the  practitioner  interested  in  solving  real-world  problems.  Commercial  im¬ 
plementations  of  similar  GUB  algorithms  have  generally  failed  to  inspire  enthusiasm 
and  are  apparently  falling  into  disuse  [Kennington  [1978]].  While  the  first  commer¬ 
cial  implementation  of  a  related  ‘'pure  network  with  side  constraints"  algorithm  )s 
(to  our  knowledge)  just  now  becoming  available,  reports  of  research  implementations 
have  generally  supported  the  view  that  such  algorithms  have  limited  application. 
The  reports  of  implementations  of  similar  "generalized  network  with  side  constraint" 
algorithms  have  offered  less  promise  than  their  pure  network  counterparts. 

■After  we  re\  iew  the  related  literature,  we  provide  a  detailed  presentation  of  the 
underl_\’ing  tableaii-ba'-ed  algoritimi  which  forms  the  foundation  of  all  our  subst-quent 


work.  We  do  so  because  the  algorithm  is  not  widely  known  and  may  be  unfamiliar 
to  the  reader.  We  then  review  the  general  factorization  approach  of  Graves  and 
McBride  [1976].  We  provide  a  design  template  for  our  developmental  approach,  and 
then  present  the  dynamic  row’  factorization  algorithms  for  GUB,  pure  network  and 
generalized  network  row  structures.  Computational  results  are  presented,  and  w’e 
then  summarize  our  conclusions  and  suggest  avenues  for  further  research. 

B.  SURVEY  OF  RELATED  LITERATURE 

While  the  terms  “partitioning”  and  “factorization”  are  frequently  used  inter¬ 
changeably  in  the  literature,  we  observe  a  distinction  between  the  two  approaches. 
We  consider  partitioning  methods  to  be  those  that  are  based  on  special  structure 
in  the  original  problem  instance.  This  structure  in  the  problem  instance  need  not 
induce  special  structure  into  the  LP  tableau,  and  in  fact  the  method  need  not  be 
tableau-based.  This  is  in  contrast  to  our  view  of  factorization,  in  which  the  algo¬ 
rithm  is  based  on  special  structure  which  occurs  in  the  simplex  basis  and  thus  in  the 
basic  tableau.  We  thus  classify  Dantzig-Wolfe  [1960]  and  Benders  [1962]  Decompo¬ 
sitions  and  Rosen's  [190-!]  f)rimal  partitioning  method  as  examples  of  partitioning 
methods. 

Perhaps  the  first  examine  of  wliat  we  consider  factorization  is  the  treatment  of 
simple  upper  bounds  by  Uantzig  [l951j  and.  independently,  by  Charnes  and  Lemke 
[1954^.  1  hey  obser\'e  that  it  is  more  efficient  to  enforce  the  “structural  "  simple 

upj)er  bound  constraints  with  logical  tests  within  the  algorithm  rather  than  treat 
them  e.xpliritl\’  along  with  other  constraints.  While  not  originally  presented  in  the 
context  of  a  formal  taiih  an  factorization,  the  approach  is  easily  viewed  as  such  and 
i:-  consistent  with  the  ucncial  aiijiroach. 
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The  mutual  primal-dual  method  of  Graves  [1965]  focuses  attention  on  the 
special  role  of  nonnegativity  constraints  in  linear  programming.  A  cle^ir  distinction  is 
drawn  between  the  computational  convenience  of  treating  nonnegativity  constraints 
implicitly  rather  than  explicitly  and  the  unambiguous  mathematical  equivalence 
of  all  problem  constraints,  structural  or  nonnegativity.  Emphasizing  the  special 
importance  of  inequality  constraints,  the  approach  yields  an  elegant  theory  (see 
Graves  [1987])  and,  as  we  wdll  see,  efficient  implementations.  We  view  this  algorithm 
as  the  first  formal  example  of  factorization. 

Dantzig  and  Van  Slyke  [1967]  extend  the  factorization  approach  applied  earlier 
to  simple  upper  bounds  in  a  more  structured  manner  in  their  treatment  of  general¬ 
ized  upper  bounds  (GUB).  In  a  problem  with  p  GUB  constraints  and  m  structural 
constraints,  their  approach  requires  a  working  basis  inverse  of  dimension  (m  +  1),  a 
considerable  savings  when  p  is  large. 

Hartman  and  Lasdon  [1972]  specialized  this  approach  to  the  multicommodity 
capacitated  transshipment  problem.  In  this  case,  the  structure  of  the  basic  pure 
network  columns  introduces  additional  structure  into  the  working  basis,  allowing 
further  simplifications  in  basis  representation  and  update  techniques.  Helgeison  and 
Kennington  [1977]  develop  techniques  for  representing  the  working  basis  inverse  in 
product  form  and  provide  graphic  interpretation  of  the  graph  updates.  Kennington 
[1977]  reports  an  implementation  of  the  algorithm. 

McBride  [1972j  and  Graves  and  McBride  [1976]  formalize  and  generalize  the 
factorization  approach.  1  hey  view  it  is  as  a  unifying  framework  for  tableau-based 
simplex  specializations  and  illustrate  this  by  developing  a  variation  of  the  GUB 
algorithm  of  Dantzic  and  \'an  Slyke  [1967]  and  a  GUB  algorithm  for  doubly  coupled 
linear  programs  of  Hartman  and  Lasdon  [1970].  They  present  a  new  algorithm 
for  the  set  t)art  it  ion  inp  liiif’ar  programming  problem  and  an  equality-constrained 
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form  of  the  pure  network  with  side  constraints  model.  McBride  [1972]  reports  an 
implementation  of  a  GUB  factorization. 

Schrage  [1975]  extends  the  approach  of  simple  and  generalized  upper  bounds 
by  considering  variable  upper  bounds  (VUB),  which  are  constraints  of  the  form 
S.  Ik,  where  i*  is  said  to  be  the  variable  upper  bound  of  Xj.  His  algorithm  allows 
the  implicit  representation  of  the  VUB  constraints  by  expressing  VUB  variables  in 
terms  of  other  basic  columns.  This  permits  the  basis  representation  to  be  treated 
in  two  parts,  one  part  a  large  matrix  which  changes  infrequently  and  thus  needs 
to  be  updated  only  occasionally,  and  the  second  part  a  small  working  basis  which 
requires  regular  attention.  Thus,  computation  and  storage  savings  may  be  realized. 
Schrage  [1978]  extends  these  ideas  to  what  he  calls  generalized  variable  upper  bounds 
(GVUB)  constraints,  which  arise  frequently  in  models  involving  fixed  charges. 

Klingman  and  Russell  [1975]  sketch  a  factorization  method  for  solving  trans¬ 
portation  problems  with  side  constraints.  They  suggest  techniques  for  performing 
simplex  iterations  and  updating  the  problem  representation.  Clien  and  Saigal  [1977] 
present  a  similar  approach  for  solving  capacitated  network  flow  problems  with  ad¬ 
ditional  linear  constraints.  Both  the  above  presentations  consider  a  graph-theoretic 
view  of  the  basis  update  mechanism  and  allow  the  basis  representation  to  be  treated 
in  two  parts,  a  part  which  corresponds  to  a  rooted  spanning  tree  defined  on  the  un¬ 
derlying  graph,  and  a  general  working  basis  inverse.  Glover,  Karney.  Klingman  and 
Russell  [1978]  report  an  implementation  of  the  Klingman  and  Russell  design,  but 
one  which  (curiously)  only  accommodates  a  single  side  constraint.  McBride  [1989] 
reports  an  implementation  which  requires  the  pure  network  rows  to  be  equalities 
and  allows  more  than  one  side  constraint. 


The  problem  of  generalized  network  problems  with  side  constraints  is  ad¬ 
dressed  by  Hultz  and  Klingman  [1976].  They  present  details  for  the  simplex  price- 
out,  column  generation  and  basis  update.  Hultz  and  Klingman  [1978]  report  an 
implementation  that  (curiously)  solves  the  “singularly  constrained”  generalized  net¬ 
work  problem.  McBride  [1989]  reports  an  implementation  that  is  not  restricted  to 
a  single  side  constraint. 

The  factorization  approach  has  been  extended  by  the  consideration  of  em¬ 
bedded  structures.  Glover  and  Klingman  [1981]  consider  a  linear  program  which 
contains  embedded  pure  network  structure,  i.e.,  the  pure  network  structure  appears 
in  only  a  subset  of  the  rows  and  columns  of  the  technological  coefficient  matrix  of 
the  problem.  Their  approach  yields  an  algorithm  similar  in  spirit  to  the  algorithms 
for  the  pure  network  with  side  constraint  model,  but  the  presence  of  the  “side  vari¬ 
ables”  significantly  complicates  the  basis  representation  and  update.  They  report 
an  implementation  of  the  algorithm  but  (curiously)  restrict  the  problem  suite  to 
problems  having  no  complicating  variables. 

McBride  [1985]  treats  the  problem  of  linear  programs  with  embedded  general¬ 
ized  network  structure.  He  presents  methods  for  pricing,  column  generation,  basis 
representation  update  and  data  structures.  A  successful  implementation  is  reported 
which  is  approximately  five  times  faster  than  MINOS  [1977]  for  the  models  tested. 

Interest  in  developing  algoritlims  to  solve  problems  with  special  substructures 
has  been  accompanied  by  work  to  identify  such  substructures  in  problem  instances. 
Greenberg  and  Rarick  [1974]  and  Brown  and  Thomen  [1980]  develop  algorithms 
to  identify  GI  B  sets.  Brown  and  Wright  [1984]  develop  algorithms  for  identifying 
pure  network  constraint  substructures.  Brown.  McBride  and  Wood  [1985]  present 
a  method  for  locating  generalized  network  structures,  both  embedded  and  row-only 


(. 


structure.-. 


Todd  [1983]  examines  factorization  from  a  geometric  standpoint  aind  constructs 
a  geometric  interpretation  which  is  in  large  measure  equivalent  to  the  algebraic 
development  of  Graves  and  McBride  [3976]. 


II.  MUTUAL  PRIMAL-DUAL  METHOD 


A.  INTRODUCTION 

This  chapter  presents  the  mutual  primal-dual  linear  programming  method  in¬ 
troduced  by  Graves  [1965]  which  provides  the  algorithmic  framework  and  notational 
conventions  for  the  research  which  follows. 

We  begin  with  a  complete  algebraic  development  of  the  primal  algorithm  fol¬ 
lowed  by  a  less  detailed  symmetric  discussion  of  the  corresponding  dual  algorithm. 
We  conclude  with  a  unified  treatment  of  the  two  which  establishes  the  theoretical 
importance  of  the  algorithm  and  justifies  its  use  as  the  foundation  of  our  specializa¬ 
tions. 

The  following  presentation  scrupulously  restates  the  Graves  [1965]  algorithm, 
incorporates  later  discussion  by  McBride  [1972],  and  accommodates  large-scale  im¬ 
plementations  by  Brown  and  Graves  [1975].  The  view  presented  here  is  not  available 
in  standard  reference  texts,  and  is  included  in  the  interest  of  completeness. 

B.  PRIMAL  PROBLEM  STATEMENT  (PLP) 

The  traditional  statement  of  the  finear  programming  (LP)  problem  is: 

(LP  )  min  :  wy 
y 

s.t.  a,y  <  r,  ,  i  =  1, .  . .  ,r7? 

>  0  ,  j  =  1,. . . ,n. 

where  y  is  an  7)-vector  of  decision  variables,  w  is  an  n-vector  of  cost  coefficients, 
each  a,  is  an  7?-vector  of  technological  transformation  coefficients,  each  r,  is  a  scalar 
right-hand  side  coefficient,  and  Cj  is  the  unit  vector.  While  this  statement  of 


the  problem  is  clear  and  unambiguous,  there  are  reasons  for  preferring  an  alter¬ 
native.  The  insistence  upon  drawing  a  formal  distinction  between  the  “structural” 
constraints  a.y  <  r,  and  the  “nonnegativity”  constraints  e^y  >  0  obscures  the  math¬ 
ematical  structure  of  the  problem  by  suggesting  that  the  two  types  of  constraints 
are  inherently  different.  Certainly  the  exploitation  of  the  special  structure  of  the 
CjV  >  0  constraints  leads  to  computational  efBciencies  in  the  implementation  of  the 
algorithm.  However,  in  the  theoretical  development  of  the  algorithm,  we  prefer  to 
treat  them  simply  as  general  inequality  constraints. 

In  order  to  achieve  a  consistent  form,  we  rewrite  the  nonnegativ'ity  constraints 
~^jy  ^  0  and  group  them  with  the  structural  constraints.  The  problem  statement 
then  becomes: 

(PLP)  min  :  icy 
y 

s.t.  :  a,y  <r,  ,  i  =  1, - m  +  n, 

where  wy  is  called  the  extremal  function.  The  constraints  of  (PLP)  define  a  set  of 
feasible  solutions  which  we  shall  call  the  feasible  set.  F: 

f  =  {y  £  /?''  1  a,y  <  r,,  i  —  I, ...  .m  +  n}  . 

Since  F  is  the  intersection  of  a  finite  number  of  closed  half-spaces.  F  is  itself  a  closed 
set.  If  F  is  nonempty  and  bounded,  then  an  optimal  solution  to  (PLP)  will  occur 
at  an  extreme  point  of  F. 

.\  point  y”  ^  /■  IS  said  to  be  a  jtasibU  point  or  feasibU  solution.  If,  for 
constraint  i.  a,]/'  <  r,  .  constraint  ?  is  satisfied  at  y”  .  and  the  quantity  r,  —  a,y°  is 
the  slack  in  constraint  i  at  y'  .  If.  on  the  other  hand,  for  constraint  i.  a,y°  >  r,  , 
constraint  i  is  violated  at  the  point  .  (he  magnitude  of  the  violation  being  a, y°  —  r, 
(the  negative  of  the  slack). 

t) 


A  point  i/°  6  /T*  is  defined  to  be  a  baste  solution  of  (PLP)  if  there  exists  an 
independent  subset  . . .  ,ai„}  of  {ai,a2)  •  •  •  ^ such  that  a,^y°  =■  r,^  for 

j  =  1, . .  .n.  Such  an  independent  subset  •  •  ,ai„}  of  {01,02, . . .  ,Om+n}  is  ^ 

basis  for  i?"  ,  and  j/°  is  the  (basic)  solution  to  the  system  Oj^y®  =  r,^,  j  =  1, . . .  ,n. 

Each  basic  solution  of  (PLP)  corresponds  to  an  extreme  point  of  F  and  for  each 
extreme  point  of  F  there  is  at  least  one  corresponding  basic  solution  of  (PLP).  Since 
there  are  at  most  (  )  ways  of  choosing  an  independent  subset  of  n  vectors  from 

{oi,  02, . . . ,  Om+n},  the  number  of  basic  solutions  of  (PLP)  and  thus  the  number  of 
extreme  points  of  F  is  finite.  Hence,  (PLP)  can  be  solved  by  searching  among  its 
basic  solutions.  The  algorithm  to  be  developed  here  will  implicitly  enumerate  the 
set  of  basic  solutions  of  (PLP)  and  terminate  in  one  of  three  states: 

1.  F  =  0 

(no  feasible  solution  exists)  : 

2.  the  extremum  is  unbounded 

(for  every  real  number  a  there  is  a  point  y°  G  F  with  u-'y°  <  o  );  or 

3.  there  exists  at  least  one  optimal  solution 

(a  point  y*  G  F  with  try"  <  u’y  V  y  G  F  ). 

We  will  first  consider  the  problem  of  finding  a  basic  feasible  solution  to  (PLP). 
Having  achieved  this,  we  will  then  consider  the  task  of  finding  a  basic  feasible  solution 
which  is  also  optimal. 

C.  OBTAINING  FEASIBILITY 

Since  we  have  included  the  nonnegativity  constraints  —e^y  <  0.  y  =  1 . n 

in  our  structural  constraint.>^  o,y  <  r,  .  7  ~  1 . 77?  -f  n,  and  since  the  origin  is  the 

unique  solution  to  the  independent  system; 
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-Cjy  =  0,  j  = 

the  origin  is  a  basic  solution  that  is  always  immediately  available  for  (PLP). 

Let  be  any  basic  solution  to  (PLP).  By  definition,  there  are  at  least  n  linearly 
independent  constraints  which  are  exactly  binding  at  y°  .  Among  the  m  remaining 
constraints,  typically  some  wiU  be  satisfied  at  y°  and  others  will  be  violated.  Our 
strategy  will  be  to  focus  on  one  violated  constraint  at  a  time,  which  we  will  call 
the  target  constraint.  Moving  from  one  basic  solution  to  another,  we  will  attempt 
to  reduce  the  violation  of  the  target  constraint  until  it  becomes  satisfied.  We  will 
restrict  our  choices  of  basic  solutions  in  such  a  way  that  all  constraints  that  are 
already  satisfied  at  y°  remain  satisfied  at  each  subsequent  basic  solution.  Once  the 
target  constraint  becomes  satisfied,  we  then  select  some  other  violated  constraint 
as  the  new  target  constraint,  and  repeat  the  process.  We  proceed  until  either  all 
constraints  are  satisfied  and  we  have  obtained  a  basic  feasible  solution,  or  we  find 
a  \iolated  constraint  which  cannot  be  satisfied,  in  which  case  we  conclude  that  no 
basic  feasible  solution  exists. 

To  formalize  these  ideas,  define  5(y°)  to  be  the  set  of  indices  of  all  constraints 
satisfied  at  a  basic  solution  : 


5(j/°)  =  |l  <  ?  <  m  +  n  I  a,y°  <  r,|  . 

Of  course.  |  S{y^’ )  |  >  n. 

Suppose  constraint  k  is  violated  at  the  basic  solution  y'^ .  Then  QkV^  >  r^  and 
k  ^  S(y^ I-  .A  necessary  condition  for  the  existence  of  a  basic  feasible  solution  is  that 
there  exists  a  basic  solution  y  with: 
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S{y)  D  S{y°) 


(2.1) 


and 

akV  <  (2.2) 

If  there  exists  no  such  basic  solution  y  satisfying  (Eq.  2.1)  and  (Eq.  2.2  )  we  conclude 
that  F  =  0.  Thus,  we  may  restrict  our  attention  to  basic  solutions  satisfying  (2.1) 
and  (2.2). 

Let  {a,,  . . .  ,a,„}  be  a  basis  for  i?"  at  j/°.  For  notational  convenience  we 

will  partition  the  constraints  into  two  sets,  those  that  are  basic  at  y°  and  those  that 
are  nonbasic  at  y°: 


Define  P  and  V  to  be  the  set  of  indices  of  the  rows  of  B  (basic  constraints) 
and  D  (nonbasic  constraints)  respectively.  Using  this  notation,  the  current  basic 
solution  may  be  ex[)ressed  as  By^  =  /,  and  since  the  rows  of  B  are  by  definition 
linearly  independent.  /?“'  exists  and  =  B~^f. 

Let  y  be  any  other  point  of  B'^‘ .  Then  y  can  be  e.xpressed  as; 


We  have  chosen  this  representation  since  at  y°  the  basic  constraints  are  satisfied 
exactly  and  thus  it  is  the  direction  vector  (y  —  y°)  of  any  proposed  move  from  y° 
to  y  that  determines  whether  or  not  a  given  basic  constraint  will  remain  satisfied 
at  y.  (  Of  course,  it  is  the  magnitude  of  (y  —  y°)  that  is  important  in  determining 
whether  the  nonbasic  constraints  remain  satisfied  as  well.) 

Now,  for  a  given  arbitrary  basis  . . .  ,p"}  of  and  any  point  y,  there 

must  exist  scalars  Aj.  . , . ,  A„  such  that; 


y  =  +  (y  -  y°)  =  y°  +  ^  A_,y  =  y°  +  PA  . 

j=i 

Since  the  choice  of  P  is  arbitrary,  let  us  choose  P  =  —B~^,  which  we  call  the 
conjugate  row  basif^.  Then: 


By  =  P(,v‘’  +  PA)  =  P(i/’  -  P-'A)  =  PyO  -  A  =  /  -  A. 

and  thus  y  satisfies  the  basic  constraints  By  <  /  if  and  only  if  A  >  0. 

Thus,  the  set  of  points  in  P’’'  which  satisfy  the  basic  constraints  b^y  <  /,  for 
?  €  B  can  be  characterized  as: 


V  e  l\''  '  V  =  V 


jei’ 


A,  >  0,;  e  B 


with  P  =  —B  ‘  . 

It  then  follows  that  every  point  y  satisfying 


can  be  expressed  as 


y  =  y"  +  PA, 

for  some  A  >  0  and  P  = 

To  address  condition  (2.2),  assume  that  we  have  selected  constraint  k  as  our 
target  constraint  (and  thus  k  is  violated  at  ,  i.e.,  >  y*,)-  To  reduce  the 

violation  of  the  target  constraint,  we  seek  a  point  y  such  that; 


dkV  <  (2.6) 

holds.  Since  such  a  point  must  also  satisfy  (2.1).  it  must  have  a  representation  of 
the  form  y  =  y°  +  PA  for  some  A  >  0.  Replacing  y  by  j/°  +  PA  .  we  have; 

PA,  <  dky\ 

wliich  holds  if  and  ordy  if  <  0.  Since  A  >  0.  it  follows  that  a  necessary 

condition  for  the  existence  of  a  point  y  satisfying  (2.5)  and  (2.6)  is  that  at  least 
one  element  of  the  vector  d,,. P  be  negative,  i.e..  <  0  for  some  J.  1  <  .?  <  n.  If 

dkP  ^  0.  we  conclude  that  F  =  0. 

Sn[)pose  that  d;  j>'  <  (i.  1  hen  (2.6)  holds  at  all  [)oints  of  the  form: 

y  =  y*’  +  A;/)'.  A,  >  (i. 

The  generic  point  i/*'  -j-  A, 7/  lies  along  an  edge  of  the  set  of  all  points  satisfying  the 
constraint'-  which  are  basic  at  y^  .  If  there  is  more  than  one  I  for  which  d^p‘  <  0. 
any  one  car,  be  cho.'-eri.  for  the  purpose  of  ex[)osition.  we  will  desicnate  the  first 
■-urb  ind'-:-:  enriiuniorrd  a--  ’ 

1  ) 


The  violation  in  constraint  k  decreases  linearly  as  A;  increases.  Constraint  k 
becomes  satisfied  when: 


dfc(y”  +  A,p')  =  gk, 

or  when  A^  =  t,  with 

A  {gk  -  dky°) 

dkP^ 

A  geometric  interpretation  is  that  y  =  +  tp‘  is  the  point  at  which  the  ray 

{y  €  I  y  =  +  A^p'.  A,  >  oj  pierces  the  hyperplane  {y  e  /?*'  |  dky  =  gk]- 

Define  A*  to  be  our  ultimate  choice  for  A;  .  Choosing  i  as  the  value  for  A* 
satisfies  the  target  constraint  k.  but  we  are  required  by  (2.5)  to  continue  to  satisfy 
all  nonbasic  constraints  which  are  alread>  satisfied  at  y°  (if  any).  The  choice  of  / 
may  cause  violations  of  such  constraints.  Thus,  Af  must  also  satisfy  the  condition: 


dj ( y"  +  A;/)' )  <  y ,  V  j  eVn  S(y^). 


Writing  this  as 


Xid.p'  5  <7;  —  dji/  . 

we  fiiul  that  we  inii'-t  choose  A;  <  ,s.  where 


lUl'i 

;Cr'.s(!.’  ) 


(y,-  -  d,y") 

"'j/' 


dj]>  >  0 


If  this  set  is  empty.  (leii:i'  Of  course,  it  is  possible  for  s  to  be  eqtial  to  zero. 
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A/  =  min{/,s},  (2.9) 

where  t  and  s  are  as  defined  in  (2.7)  and  (2.8)  respectively,  we  obtain  the  largest 
reduction  in  the  violation  of  the  target  constraint  consistent  with  the  fe^lsibility 
restriction  (2.5)  with  respect  to  the  chosen  direction  p'  . 

The  selection  A*  according  to  (2.9)  leads  to  a  new  basic  solution  =  y°  +  A*p*,. 
If  A*  >  0  at  each  step,  the  transitivity  of  (2.5)  and  (2.6)  guarantees  that  no  basic 
solution  will  be  repeated  before  any  given  target  constraint  is  satisfied,  or  until 
the  conclusion  is  reached  that  F  =  0.  At  any  given  step,  either  the  'easible  set  is 
shown  to  be  empty,  a  basic  solution  is  found  which  satisfies  the  target  constraint 
or  a  positive  reduction  in  the  violation  of  the  target  constraint  is  achieved  as  the 
result  of  moving  to  another  basic  solution.  Since  the  set  of  basic  solutions  is  finite, 
the  third  alternative  may  be  repeated  a  finite  number  of  times  before  one  of  the 
first  two  alternatives  occurs.  If  this  constraint  becomes  satisfied,  (2.5)  guarantees 
that  every  subsequent  solution  will  also  satisfy  the  constraint.  We  may  then  select 
a  new  violated  nonbasic  constraint  as  the  target  constraint.  Since  there  are  a  finite 
number  of  constraints,  we  will  either  discover  a  basic  feasible  solution  or  determine 
in  a  finite  number  of  steps  that  none  exists. 

If  \]  =  0  at  any  step  of  the  algorithm,  we  say  that  a  block  has  occurred. 
Blocks  will  be  discussed  later  in  detail. 

D.  FROM  FEASIBILITY  TO  OPTIMALITY 

Once  feasibility  is  achieved,  say  at  .  the  process  of  proceeding  to  optimality 
can  be  thought  of  as  minimizing  over  F  the  violation  of  the  constraint  v:y  <  wy^  —  a. 
where  o  is  a  sufhciently  large  [xisitive  constant. 

It; 
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Let  t/°  be  a  basic  feasible  solution  with  By°  =  f  and  P  =  —B~^  .  Since  every 
point  y  G  i?"  which  satisfies  By  <  f  may  be  written  as  y  =  y°  +  P\  for  some  A  >  0, 
the  value  of  the  extremal  function  wy  at  such  a  point  y  is: 

n 

rvy  =  w{y°  +  PA)  =  wy°  +  wPX  =  wy^  +  \jWjP . 

i=i 

Thus,  a  necessary  condition  for  the  existence  of  a  point  y  £  F  such  that  xvy  <  wy^ 
is: 

u'p-’  <  0  for  some  j.  1  <  j  <  n. 

If  wfP  >  0,  then  y°  must  be  a  feasible,  optimal  solution  to  (PLP). 

Suppose  U'p'  <  0.  Since  all  constraints  are  satisfied  at  y°,  the  greatest  reduc¬ 
tion  in  the  value  of  the  extremal  function  in  the  direction  is  achieved  when  A*  is 
chosen  to  be  s.  where  s  is  as  defined  in  (2.8).  If  5  =  oc  ,  (PLP)  has  an  unbounded 
extremum.  If  0  <  Af  <  oc  .  a  positive  reduction  in  the  value  of  the  extremal  function 
can  be  achieved  by  moving  from  the  basic  feasible  solution  y°  to  the  basic  feasible 
solution  y  =  y°  4-  Afp^  .  If  X]  =  0.  a  block  is  encountered. 

Once  a  feasible  solution  has  been  obtained.  (2.5)  and  (2.6)  become  equivalent 
to: 

y  G  P  .  (2.10) 

try  <  icy°  .  (2.11) 

The  transitivity  of  (2.10j  and  (2.1 1  )  ensures  that  no  basic  feasible  solution  will 
ie-  repeated  a:^  lone  a-  A"  >  0  at  each  step.  .At  each  iteration,  either  a  basic  feasible 
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solution  is  found  to  be  optimal,  the  solution  is  found  to  be  unbounded  or  a  positive 
reduction  in  the  value  of  the  extremal  function  is  achieved  by  moving  to  another 
basic  feasible  solution.  Since  the  set  of  basic  feasible  solutions  is  finite  and  no  basic 
fetisible  solution  is  repeated,  the  third  alternative  can  only  occur  a  finite  number  of 
times.  Thus,  if  A*  >  0  at  each  step  and  if  a  finite  optimal  solution  exists,  it  will  be 
found  in  a  finite  number  of  steps. 

E.  BLOCKS 

Suppose  that  y°  is  a  basic  solution  satisfying 


B/  =  /. 

where,  as  before,  B  is  of  dimension  n  by  n  and  nonsingular  and  is  the  column 
of  P  =  — P“*.  Then.  y°  is  the  unique  point  in  P"  lying  in  the  intersection  of  the  n 
hyperplanes; 


€  R''  I  =  /,.  ?  =  1 . Tij. 

Suppose  that  p’  is  a  direction  that  either  leads  to  a  reduction  in  the  violation 
of  some  target  constraint  or,  if  ?/“  is  feasible,  a  reduction  in  the  value  of  the  extremal 
function.  Since  the  nonbasic  constraints  arc  of  the  form: 


J  €  P. 

a  block  occurs  whc!;.  fo:  at  least  one  nonbasic  constraint 


=  p,.. 


and 


dkp‘  >0.  (2.12) 

Thus,  any  movement  away  from  y°  in  the  direction  p‘  will  result  in  violation  of  the 
nonbasic  constraint  dky°  <  gk  ■ 

If  d^y^  =  gk  ,  then  is  in  a  sense  “over-determined”,  and  i/°  is  said  to  be 
a  “degenerate”  solution.  Geometrically,  t/°  lies  in  the  intersection  of  at  least  n  -|-  1 
hyperplanes.  Algebraically,  there  is  more  than  one  basis  that  can  be  formed  from 
the  row  vectors  {ci,  qj,  . . . ,  Um+n}  for  which: 


holds,  and  there  is  more  than  one  basic  solution  which  corresponds  to  the  extreme 
point  j/°  .  A  block  is  therefore  encountered  when  an  over-determined  solution  sat¬ 
isfies  (2.12)  in  an  “improving”  direction  (e.g.,  one  that  leads  either  to  a  reduction 
in  infeasibility  or.  if  feasible,  to  a  reduction  in  the  value  of  the  extremal  function). 
W’e  will  develop  a  method  for  dealing  with  blocks  later. 

F.  THE  PRIMAL  ALGORITHM  AND  A  SUPPORTING  BASIC  TAB¬ 
LEAU 

The  primal  algorithm  proceeds  as  follows: 

1.  Identify  an  initial  basic  solution.  Notice  that  the  origin  satisfies 


and  thus  may  always  serve  as  the  initial  solution. 


lb 


2.  If  the  current  solution  is  infeasible,  select  a  violated  primal  constraint  index  k 

(i.e.,  an  index  k  for  which  <  0).  This  requires  the  quantity: 

g-Dy^  . 

3.  Within  the  target  row,  select  an  element  of  the  proper  sign  (negative,  according 
to  our  convention)  whose  index,  /,  specifies  a  “traxisformation  column”.  If  the 
current  solution  is  infeasible,  this  requires  the  quantity: 

dkP  , 

with  the  transformation  column  satisfying  d^p^  <  0.  If  the  current  solution  is 
feasible,  this  requires  the  quantity: 

wP 

with  the  transformation  column  satisfying  u>p'  <  0.  If  no  such  element  exists, 
then  the  problem  is  infeasible  (if  the  current  solution  is  infeasible  and  di^P  >  0) 
or  optimal  (if  the  current  solution  is  feasible  and  wP  >  0).  The  task  of 
selecting  such  an  index  I  is  commonly  referred  to  as  “pricing"  or  a  “pricing 
strategy" . 

4.  Compute  t  as  in  (Equation  2.7).  If  the  current  solution  is  feasible,  assign 
i  =  oc. 

5.  Compute  s  as  in  (Equation  2.8).  This  computation  is  commonly  called  a  “ratio 
test",  or  a  “minimum  ratio  test".  This  computation  requires  the  quantities: 

g  —  Dy^  and  Dp‘ . 

G.  Compute  A:  a^  in  (Equation  2.9).  If  At  =  oc.  the  problem  is  unbounded. 


7.  Update  the  current  solution  according  to  the  computation; 

y’  =  t/°  +  Arp'. 

8.  Update  the  assignments  of  constraint  indices  to  B  and  T>,  and  update  P. 

9.  Go  to  Step  2. 

The  only  quantities  needed  at  each  step  of  the  algorithm  are  the  matrix  DP^ 
the  column  vector  g  —  Dy^  and  the  row  vector  wP.  For  notational  convenience,  we 
can  define  D  to  be  an  (m  +  1)  by  n  matrix  whose  bottom  row  is  w,  and  g  to 

be  a  (m  +  1)  dimension  column  vector  whose  bottom  entry,  Pm+i,  is  zero.  Then 

yrr.+i  -  dm+iy'’  =  0  -  wy°  =  -U'y°, 

which  is  the  negative  of  the  extremal  function  value  at  the  point  y°  .  We  can  now 
conveniently  display  all  the  relevant  information  by  forming  the  (m  +  1)  by  (m  +  1) 
matrix: 

[DP\  9 -by].  (2.13) 

which  we  will  call  the  basic  tableau. 

By  computing  the  basic  tableau  in  partitioned  matrix  form,  we  may  isolate  the 
important  algcl.>raic  components  required  by  this  method.  To  do  so.  let  us  assume 
that  at  the  current  ba^-ic  solution  the  basis  consists  of  h  structural  constraints  and 
(7i  —  h)  nonnegat i vit y  constraints.  Then  (2.3)  can  be  written  as; 


■  b, 

a.i 

■/i 

’  '•n  ' 

62 

/2 

bn 

bh+i 

= 

a.* 

and  /  = 

A 

fh-^1 

= 

0 

bh+2 

“®Jh+2 

fh+2 

0 

. 

.  A 

0 

In  partitioned  matrix  form, 


j 

and  thus 


\  0 


n— /i 


Aj2  J  h 
-If  }  n  -  /i  , 


P  = 


h 


\  0 


n-/i 


/  J 


}h 

)n 


h 


Similarly.  (2.4)  can  be  written  as: 


■  d,  ■ 

-tji 

■  0 

fl'2 

-^2 

92 

0 

J 

and  g  = 

9h 

0 

9h+l 

^n.+i 

dh+i 

^>h*2 

9h+2 

^•h+2 

d.:. 

^1,,  +  r,, 

.  9rr. 

*n-+  Til 

22 


(2.14) 


and  in  partitioned  matrix  form: 


h 


n— A 


-I  0  }h 

^21  ^22  }m  —  h 


(2.15) 


Then 


DP  =  -DB-^  = 


atM 


11  ^12 


\~-^2i>4j/  A22  —  .^2i-^ii''^i2  /  }  rn  —  h  , 

and  we  shall  call  DP  the  principal  part  of  the  tableau.  By  partitioning  w  = 
(u;i,u;2),  =  (51,^2)^  and  =  iVuV^)^  basic  tableau  (2.13)  may  be  written 

in  partitioned  matrix  form  as: 


h 

n  —  h 

1 

Au' 

All  Ai2 

y? 

)h 

-A21  Aj~]' 

A22  —  A2iAj]Mi2 

.<72  -  A2iyi  - 

^22^2 

}  m  —  h 

-u-iAfi' 

u'2  -  u’lAn  Ai2 

0 

1 

}l  • 

(2.16) 


Note  that  j/°  is  displayed  explicitly  in  the  tableau  (2.16).  Also,  ^2  =  0  since  the 
corresponding  nonnegatii  ity  constraints  are  basic  and  thus  binding. 

As  the  algorithm  proceeds  from  one  basic  solution  to  another.  (2.16)  can  be 
updated  using  a  slight  variant  of  the  Gauss-Jordan  transformation  (which  we  will 
call  a  pivot).  The  transformation  is  as  follows: 

Let  \bP\g  —  by°]  =  [r',j]  be  the  basic  tableau  at  any  basic  solution  .  If 
=  dsp'  ^  then  the  basic  tableau  al  the  basic  solution: 


1  ('  .  /  I  \  I  (I 


(g.-:  -  I 


written  as 


D'P'  I  g'  -  D'y']  = 


is  calculated  as: 


general  elements  {i  ^  s^j  ^  1) 

pivotal  column  (i  s,  j  =  1)  •. 

,  _  —v,i 

pivotal  row  {i  =  s,  j  1)  : 

*’•'  VsJ 

pivotal  element  t  =  s,  j  =  I 
,  _  1 

After  applying  the  transformation  to  (2.16),  row  and  column  interchanges  may 
be  necessary  to  restore  the  explicit  structure  of  B  in  (2.14)  and  D  in  (2.15).  A  proof 
of  this  may  be  found  in  Graves  [1987].  who  also  extends  the  basic  algorithm  to 
include  free  variables,  equality  constraints  and  bounded  variables. 


G.  DUAL  LINEAR  PROBLEM  (DLP) 

Corresponding  to  (PLP)  is  the  following  linear  program  called  its  dual: 

(DLP)  m.ax  :  xr 

T 

s.t.  :  xA  <  w 
xl  <  0. 

The  relationship  betwf'<ni  (DLP)  and  (PLP)  is  as  follows: 
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1.  if  (PLP)  is  infeasible,  then  (DLP)  is  either  infeasible  or  unbounded. 

2.  if  (PLP)  is  unbounded,  then  (DLP)  is  infeasible. 

3.  if  (PLP)  possesses  a  finite  optimal  solution  y’  ,  then  (DLP)  also  possesses  a 
finite  optimal  solution  i*  and  wy“  —  x'r. 

See  Graves  [1987]  for  a  proof  of  this. 

Solving  (DLP)  directly  has  important  advantages  in  certain  problems,  and  we 
will  show  that  its  solution  plays  an  important  role  in  dealing  with  blocks  in  (PLP). 
We  will  thus  develop  a  dual  algorithm  in  a  manner  similar  to  that  for  (PLP),  sparing 
some  of  the  symmetrical  detail  where  possible. 

Let  Xo  be  a  basic  solution  for  (DLP).  Then  there  exists  an  independent  subset 

{a-” , a-’%  . . . , }  of  _ such  that  Xoa^'  =  w^\  i  =  I, ...  ,m,  and 

this  linearly  independent  subset  will  be  called  the  basis  for  R'^  at  Xo  . 

Partitioning  the  dual  constraints  into  those  which  are  basic  at  Zo  and  those 
which  are  nonbasic  at  yields: 

T  =  {i\t\...,r)  r=  (a^^a^^...,a''") 

u  =  {v\v^ . ir)  =  (U’^^U•^^...,U’'"’),  (2.17) 

and 

K  -  . H-'" )  =  (a-’"’+‘ ,  . . . ,  0-’'"+" ) 

n  =  (z'.r^ . f")  =  (»•■'">+’,  (2. IS) 

Define  T  and  K  to  be  the  set  of  indices  of  the  columns  of  T  (basic  dual 
constraints)  and  A  inwnbasic  dual  constraints),  respectively.  The  current  basic 
solutic.’r:  x._  nia\'  tli'M;  be  e.xpre'^sed  a'- 


Xo  —  uT  \ 


and  any  other  point  i  in  may  be  expressed  as 

I  =  +  (x  -  Xo). 

For  a  given  arbitrary  basis  {91,92?  •  •  •  ?9m}  of  ,  any  such  point  x  may  be  repre¬ 
sented  as: 

m 

^  +  =  lo  +  rl'Q- 

t=i 

The  current  basic  constraints  are  XoO^'  -  w^' and  for  these  con¬ 
straints  to  remain  satisfied  at  the  generic  point  x  =  x^  -t-  %'Q.  we  must  have 

[la-’’  <  =  1 . m].  Choosing  Q  =  T~^  as  a  convenient  basis  for  at  x^.  we 

then  have: 

xT  =  (xo  +  n:Qyr  =  xoT  -f  pQT  =  u  -f  C, 
and  thus  we  must  have  v  <  (1. 

Now  consider  a  violated  dual  nonbasic  constraint  t.  where  a  violation  means 

x,C  >  r'. 

To  reduce  the  violaiion.  we  seek  a  point  x  at  which 

xA’'  <  X.  k' . 


holds.  Such  a  point  must  also  satisfy  the  basic  constraints,  and  thus  the  condition 
for  a  reduction  in  the  violation  of  the  dual  target  constraint  is: 

xk*  =  (xo  +  rl^Q)k*  =  Xok*  +  ipQk^  <  Xok\ 

which  implies  that 

i^Qk*  <  0. 

Since  V  <  0.  a  necessary  condition  for  constraint  t  to  be  satisfied  is  that  there  exists 
an  i  with  q^k^  >  0.  If  none  of  the  elements  of  q,k^  are  positive,  we  conclude  that  the 
dual  constraints  are  inconsistent. 

The  dual  algorithm  proceeds  along  the  edges 

X  =  Xc,  +  v’q^. 

from  basic  solution  to  basic  solution.  We  define 

=  |l  <  j  <  ni  -r  n  \  <  u-’l  . 

and  insist  that 

r)  D  5(>J. 

To  satisfy  thi'^  condition,  we  consider  the  effect  of  moving  along  the  edge  of  a 
general  nonbas'c  dual  coii'^traint  I:  which  i:"  currently  satisfied  at  .  We  must  have 


xk^  =  {xo  +  xl>'g,)lr’  =  Xgk^  +  ^  ^■’• 


Since  the  constraint  is  satisfied  at  Xo  ,  we  have 

Xok^  <  tr'  =>  tr’  —  Xok^  >  0, 

and  since  <  0,  if  g,k^  <  0  then  as  0’  decreases  in  value  we  approach  the  boundary 

of  the  constraint.  If  —  Xok^  =  0  then  <  0  causes  an  immediate  violation  if 

tp'  <  0.  Thus,  a  block  has  occurred. 

We  have  shown  that  our  choice  for  t/’*,  denoted  v*.  should  be 

V'l  =  max  {6, c},  (2.19) 

where 

6  =  (r' -  (2.20) 

and 

c  =  max  {(r’  -  J,.F)/g.F  |  m’  -  x.k^  >  0.  g.F  <  oj  .  (2.21 ) 

Once  dual  feasibility  has  been  achieved,  we  wish  to  maximize: 

xj'  =  x,.r  +  f(Qr]. 

Since  f  <  0.  a  necessa’y  condition  to  achieve  an  increase  in  tlie  value  of  the  extremal 
function  is  that  there  exists  an  ?  such  that  g,r  <  0.  Hence,  when  q,r  >  0.  we  conclude 
that  x^,  is  dual  ojit  imai, 

7  he  dual  alccriiiii;;  thei,  pioceed-  as  follows; 


1.  Identify  an  initial  basic  solution.  Notice  that  the  origin  satisfies 


x°  •  /  =  0, 

and  thus  may  always  serve  as  the  initial  solution. 

2.  If  the  current  solution  is  infeasible,  select  a  violated  dual  constraint  index  t 
(i.e.,  an  index  t  for  which  r‘  —  Xok*  <  0).  This  requires  the  quantity: 

V  —  XqK. 

Column  t  is  then  referred  to  as  the  ‘‘target  column’".  If  the  current  solution  is 
feasible,  the  right-hand  side  column  is  designated  the  “target  column’". 

3.  Within  the  target  column,  select  an  element  of  the  proper  sign  (positi%'e,  ac¬ 
cording  to  our  convention)  whose  index,  i.  specifies  a  “transformation  row". 
If  the  current  solution  is  infeasible,  this  requires  the  quantity: 

Qk'  . 

with  the  transformation  row  satisfying  q^k^  >  0.  If  the  current  solution  is 
feasible,  this  requires  the  quantity: 

Qr  . 

with  the  transformation  row  satisfying  q,r  >  0.  If  no  such  element  exists,  then 
the  problem  is  infeasible  (if  the  current  solution  is  infeasible  and  q,k^  <  0)  or 
optimal  (if  the  current  solution  is  feasible  and  q,r  <  0).  The  task  of  selecting 
such  aii  index  ?  is  commonly  referred  to  as  “pricing",  or  a  "pricing  strategy"  . 

•S.  Compute  k  as  in  (Equation  2.20’!.  If  the  current  solution  is  feasible,  assign 
h  —  —  oc  . 

•>ij 


5.  Compute  c  as  in  (Equation  2.21).  This  computation  is  commonly  called  a 
“ratio  test”.  This  computation  requires  the  quantities: 

V  —  xqK  and  Qk\ 

6.  Compute  xpl  as  in  (Equation  2.19).  If  xpl  =  oo,  the  problem  is  unbounded. 

7.  Update  the  current  solution  according  to  the  computation: 

Xi  =  Xo  + 

8.  Update  the  assignments  of  constraint  indices  to  T  and  AC,  and  update  Q. 

9.  Go  to  Step  2. 

This  development  shows  that  at  each  step,  the  quantities  Q/C,  (r  —  XoK)  and 
Qr  are  needed  to  proceed  with  the  dual  algorithm. 

To  develop  a  matrix  partitioned  form  of  the  dual  tableau,  we  proceed  as  before. 
Assuming  the  dual  basis  consists  of  h  structural  constraints  and  m  —  h  nonnegativity 
constraints,  we  have: 


.,r-j  =  (o-'’ ,  a-'A  . . 

^  Am 

-  .  . 

=(u-.u-^A. 

. .  .u-"\0.0. 

....0). 

u  =  =  (u''.0). 


The  nonl)asi(  constraints  are  then: 


A  =  I . V' )  =  {f'- . f''*.  . 0'^" 


?■=!'■  .  1 


l:'‘  1  =  !f).n . 0.  !/•’'■*  A 

.'•ifl 


so  V  =  =  (0,ll>^)  . 

The  matrix  partitioned  form  of  the  basis  is  then 


I— A 


T  = 


I  \ 

Au  0 


\  A-ii 

and  with  the  choice  oi  Q  =  T~^ 
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}h  , 

}m-h  , 


tn—h 
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Q  = 


^11 
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7  J 

with  the  remaining  constraints  forming 


]h  , 

}  m  —  h  , 


/  >1 
0  A22 


The  principal  part  of  the  dual  tableau  is: 


QK 


0  ■ 

[  7  Ai2  1 

-.42, /ir,’  / 

0 

A{,’A:2 

—  .421.4i/  A22  —  A2\A^^  A-i2 


which  we  find  to  be  exactly  the  principal  part  of  the  primal  tableau,  so 


DP  =  QK  . 


The  quantitit'>  (r  —  A' i  arc 


r  -  .7  A’ 


f  -  iuQ)K 


V  —  uDP 


=  (0,u;')  q"'  All  ^12  _^j^p 

=  vP  —  uDP 
=  {v-uD)P 

=  {(o--’) Lll'’ 

=  {{0,w^)  -  (-w\0)}  P 
=  wP. 

The  quantity  Qr  is: 


=  (?{P  +  A'/} 


—  Qg  +  Qh  f 
_  ■  O' 

— '421-4  j/  1 


+  DPf 


=  g  +  D{Pf) 


0 

r2 


+  DPf 


=  g-Df. 


We  thus  find  that  the  quantities  required  for  the  dual  tableau  [QK.  r  — x^A'  and 
Qr)  are  exactly  the  quantities  required  for  the  primal  tableau  [DP.  u'P  and  g  —  Dy^). 
Therefore,  the  primal  and  dual  algorithms  use  the  same  tableau.  Reviewing  the 
summa’'ie.‘-  of  the  primal  and  dual  algorithms,  we  see  the  strong  symmetry  between 
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the  two.  Operating  on  a  single  tableau,  the  primaJ  algorithm  identifies  a  target 
row,  selects  a  transformation  column,  performs  a  ratio  test,  and  performs  a  pivot 
update,  while  the  dual  algorithm  identifies  a  target  column,  selects  a  transformation 
row,  performs  a  ratio  test  and  performs  a  pivot  update.  This  remarkable  symmetry 
allows  us  not  only  to  perform  either  algorithm  on  a  single  tableau,  but  to  switch 
between  one  edgorithm  and  the  other  as  often  as  we  choose.  This  fact  has  important 
implications  for  our  implementation  and  it  enables  us  to  deal  with  blocking  in  a 
systematic  and  consistent  way. 

H.  RESOLUTION  OF  BLOCKING 

The  resolution  of  blocking  in  either  the  primal  or  dual  algorithm  will  be  accom¬ 
plished  by  shifting  to  the  alternate  algorithm  when  blocking  occurs.  The  alternate 
algorithm  is  applied  to  a  subproblem  of  the  original  cind  at  worst  we  are  lead  to  a 
contracting  sequence  of  problems  to  which  we  alternately  apply  the  primal  and  dual 
algorithm.  A  strict  contraction  is  assured,  and  thus  in  at  most  a  finite  number  of 
steps  resolution  is  achieved. 

We  will  demonstrate  this  procedure  by  assuming  that  we  have  started  with 
the  primal  algorithm.  When  a  block  occurs,  assume  we  have  rearranged  the  rows 
so  that  the  zeros  in  the  right-hand  side  column  occur  contiguously. 

Let  A'l  =  7?  -I-  1  and  —  1  where  t  is  the  index  of  the  target  row  (if  the  current 
solution  is  primal  feasible,  =  ni  I .  the  extremal  function  row).  Let  be  the 
pivot  column  which  caused  the  block.  The  situation  is  shown  in  Figure  2.1. 

Define  Subproblem  ( 1 )  as  consisting  of  all  the  columns  of  the  original  problem, 
but  only  those  rows  with  zeros  in  the  right-hand  side  column,  and  as  shown  in  Figure 
2.2.  "extrema!  function"  row  L2.  the  target  row  of  the  original  primal  problem.  Row 
A.,  is  the  row  containing  tlie  blocking  element.  Now  apply  the  dual  algorithm  to 


Figure  2.1:  Primal  Block  in  Column  A:3 


ki 

^2 

Subproblem  (1).  Because  the  right-hand  side  column  of  any  potential  pivot  row  is 
zero,  the  right-hand  side  column  of  the  original  primal  problem  is  invariant  to  any 
pivot  in  Subproblem  (1).  Let  i',^  denote  the  element  in  row  i  and  column  j  of  the 
tableau.  We  proceed  with  the  dual  algorithm  on  Subproblem  (1)  until  one  of  the 
following  situations  occur  (which  may  be  after  zero  or  more  dual  pivots): 

L  pivot  in  column  kj  no  longer  reduces  the  violation  of  primal 

constrrint  a-j  ,  and  thus  the  primal  block  has  been  eliminated.  We  thus  return 
to  the  original  primal  problem  and  seek  a  new  columm  in  which  to  pivot. 

2.  ^'kt.ky  ^  b.  A  pivot  in  column  kj  no  longer  threatens  to  violate  primal  con¬ 
straint  .  We  thus  compute  new  ratios  using  columns  ki  and  ^-3  of  the  original 
primal  tableau  to  determine  if  we  can  make  a  gain  on  the  target  constraint  k^ 
using  column  (  .  . 

:n 


^3  ^1 


+ 

0 

0 

0 

. 

Figure  2.2:  Subproblem  (1) 


Figure  2.3:  Dual  Block  in  Subproblem  (1) 
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3.  We  encounter  a  block  in  the  dual  algorithm  applied  to  Subproblem  (1).  This 
occurs  when  a  column  of  Subproblem  (1)  contains  a  negative  element  in  row 
ki  and  a  zero  in  row  /i-2,  as  illustrated  in  Figure  2.3. 

Let  ks  be  the  column  in  Subproblem  (1)  causing  the  dual  block.  We  now 
define  a  further  contraction,  Subprobiem  (2),  and  switch  to  the  primal  algorithm. 
Subproblem  (2)  consists  of  all  the  rows  of  Subproblem  (1).  but  only  those  columns 
of  Subproblem  (1)  with  zeros  in  its  target  row  (^'2).  Note  that  the  current  target 
column  in  Subproblem  (1)  (^-3)  must  have  a  negative  element  in  its  “extremal"  row 
(^2)'  As  a  notational  convenience  we  reverse  the  sign  of  this  element  and  record  this 
action  by  labeling  the  column  as  — L3.  Column  —k^  becomes  the  right-hand  side 
column  of  Subprobiem  (2).  as  shown  in  Figure  2.4. 

The  number  of  columns  in  Subproblem  (2)  is  reduced  by  at  least  one  from 
the  number  of  columns  in  Subprobiem  (1).  We  now  apj)ly  the  primal  algorithm  to 
Subproblem  ^2)  until  (which  may  be  after  zero  or  more  primal  pivots); 
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Figure  2.4:  Subproblem  (2) 
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1.  Vkt-ki  >  0,  implying  inat  <  0.  A  pivot  in  column  k^  no  longer  threatens 
to  violate  primal  constraint  k4  .  We  have  thus  removed  the  original  primal 
block,  and  we  proceed  by  returning  to  the  original  primal  problem,  and  com¬ 
puting  new  ratios  using  columns  k^  and  k^  . 

2.  >  0.  A  pivot  in  row  no  longer  threatens  to  violate  dual  constraint 
fcs,  and  thus  ac  have  removed  the  dual  block  in  Subproblem  (1).  We  revert  to 
Subproblem  (1)  and  continue  with  the  dual  algorithm. 

3.  We  encounter  a  block  in  the  primal  algorithm  applied  to  Subproblem  (2).  We 
illustrate  this  situation  in  Figure  2.5. 

Let  fce  be  the  row  causing  the  primal  block  in  Subproblem  (2).  Since  the  the 
target  row  in  Subproblcm  (2)  (^-2)  is  identically  zero  for  all  but  the  right-hand  side 
column  (  —  kj),  primal  unboundedness  cannot  occur.  Also,  the  entire  k2  row  in  the 
original  problem  is  invariant  to  pivots  in  this  subproblem.  We  proceed  just  as  we  did 
w’hen  faced  with  this  situation  in  the  original  primal  problem.  The  current  target 
row  (^4)  in  the  current  problem  (SubprobJem  (2))  becomes  the  extremal  function 
row  in  a  newly  defined  contraction  (Subproblem  (3)).  All  columns  of  the  current 
problem  are  retained,  but  only  those  rows  with  zeros  in  the  right-hand  side  column 
of  the  (  —  k^)  current  problem  are  carried  forward  to  the  contraction.  Ihis  again 
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Figure  2.5:  Primal  Block  in  Subproblem  (2) 
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a  strict  reduction  in  the  nurnher  of  rows.  The  new  subproblem  is  shown  in 
Figure  2.6.  We  then  proceed  with  the  dual  algorithm,  exactly  as  before. 

The  construction  of  the  contracting  subproblems  through  the  nested  sets  of 
zeros  in  the  columns  and  rows  guarantees  a  monotonic  decrease  in  the  sizes  of  the 
higher-order  subproblems.  This  ensures  the  ultimate  resolution  of  degeneracy  and 
gives  us  a  complete,  symmetric  and  finite  algorithm  for  the  solution  of  LP  problems. 

The  blocking  resolution  scheme  given  here  is  a  constructive  algorithm  to  iden¬ 
tify  strictly  improved  solutions.  The  restricted  subproblems  ultimately  yield  a  pivot 
sequence  satisfying  a// higher-order  criteria.  Geometrically,  we  systematically  search 
the  degenerate  subspace  for  an  improved  representation.  This  is  in  sharp  contrast 
to  ad  hoc  “anti-degeneracy”  and  “anti-cycling”  schemes  which  invoke  arbitrary  sec¬ 
ondary  mechanisms  not  at  all  related  to  the  geometry  or  mathematics  of  the  problem 
at  hand,  and  consequently  admit  nuisance  degenerate  pivots  with  no  constructive 
motivation. 

I.  RELATIONSHIP  BETWEEN  PRIMAL-DUAL  ALGORITHM  AND 
SIMPLEX  METHOD 

We  now  depart  from  the  presentation  of  Graves  [1965]  to  discuss  the  relation¬ 
ship  between  the  Mutual  Primal- Dual  Method  and  the  classical  Simplex  Method. 
After  a  brief  discussion  of  the  similarities,  we  will  explain  our  reasons  for  adapting 
the  Primal-Dual  view, 
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Figure  2.6:  Subproblem  3 


The  Simplex  Method  assumes  primal  feasibility,  and  we  must  identify  a  start¬ 
ing  basic  feasible  solution.  Because  this  is  celdom  practical,  the  usual  approach  is  to 
formulate  and  solve  a  new  linear  program  closely  related  lo  the  original  which  has 
the  same  optimal  solution  (assuming  one  exists)  and  possesses  an  easily  idenl'^^ed 
basic  feasible  solution.  This  related  problem  is  derived  from  (PLP)  by  augmenting 
it  with  slack  variables,  resulting  m  tne  foiiowing: 


(APLP) 

min  :  wy  -|- 

Os 

s.t.  ;  Ay  + 

Is  =  r 

}  m 
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>  0 
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Is  >0 

}  m 

In  the  classical  Simplex  view,  a  basis  is  a  collection  of  m  linearly  independent 
columns.  Let  As  he  such  a  simplex  basis  which  corresponds  to  a  primal  feasi¬ 
ble  solution  for  (APLP).  Suppose  we  partition  the  coefficient  matrix  in  a  manner 
determined  by  the  basic  slack  variables,  yielding  (after  possible  row  and  column 
interchanges): 
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and  the  matrix  of  nonbasic  columns  is 


The  basic  variables  are  ys  =  (yi,'S2)  while  the  nonbasic  variables  are  y^B  = 
(-Siiyz)?  the  principal  part  of  the  Simplex  tableau  is 

Borrowing  from  the  perspective  of  the  Primal-Dual  algorithm,  we  may  generate 
the  Simplex  tableau  by  partitions: 
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and  the  table»^u  becomes 
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which  is  precisely  the  p^rt  of  the  Primal-Dual  tableau.  This  is  no  surprise 

when  we  note  that  the  basis  for  the  Dual  Algorithm,  T,  is  exactly  As,  and  the  Primal 
and  Dual  Algorithms  share  the  same  tableau.  We  may  interpret  the  Primal  and  Dual 
Algorithms  as  simply  two  different  perspectives  of  the  same  tableau,  wherein  the 
Primal  Algorithm  a  pivot  is  viewed  as  exchanging  primal  constraints  and  in  the 
Dual  Algorithm  a  pivot  is  viewed  as  exchanging  dual  constraints.  The  classical 
Simplex  Method  may  then  be  interpreted  as  solving  (PLP)  using  the 
Dual  Algorithm  perspective.  That  the  classical  Simplex  Method  is  naturally 
interpreted  as  a  dual  algorithm  comes  as  a  surprise  to  the  conventionally  trained. 
However,  the  consequent  mathematical  insight  is  compelling,  especially  in  light  of 
the  notational  simplification  and  apparent  underlying  role  of  Aj"/  . 


There  are  several  reason.'  for  preferring  the  Primal-Dual  Algorithm  to  the 
Simplex  Method.  From  a  computational  standpoint,  because  slack  variables  are 
carried  logicall\'  rather  than  introduced  explicitly,  we  are  able  to  clearly  identify  the 


essentiaJ  information  needed  to  execute  the  algorithm.  The  matrix  Af/  plays  a  key 
role  in  the  calculation  of  the  tableau,  and  the  entire  tableau  can  be  constructed  from 
and  original  problem  data.  Since  is  a  submatrix  of  ,  it  is  smaller  and 
requires  fewer  arithmetic  operations  to  update  than  does  . 

A  second  advantage  of  the  Primal-Dual  Algorithm  lies  in  the  flexibility  it 
offers  for  specialisation  to  particular  problem  classes  or  structures.  Indeed,  it  is  the 
special  structure  and  simplicity  of  the  nonnegativity  constraints  that  motivate  the 
development  of  the  algorithm  in  the  first  place.  It  is  frequently  the  case  that  other 
special  structures  can  be  identified  in  classes  of  (PLP).  Examples  of  such  structures 
include  simple  upper  bounds,  generalized  upper  bounds,  variable  upper  bounds, 
pure  and  generalized  network  substructures,  etc.  Such  structure  may  be  static  in 
that  its  nature  and  dimension  remains  fixed  throughout  the  solution  process,  or 
the  structure  may  be  dynamic  in  which  case  its  precise  nature  and/or  dimension 
may  vary  as  the  problem  is  solved.  Some  special  structures  may  be  more  strongly 
characterized  by  the'r  column  structure  and  others  by  their  row  structure.  The 
Primal-Dual  .Algorithm  allows  one  to  effectively  exploit  virtually  any  such  problem 
structure  in  a  natural  manner  and  greatly  simplifies  the  implementation  of  such  a 
specialization. 

When  the  linear  programming  problem  appears  as  a  subproblem  in  a  more 
sophisticated  solution  setting  (for  example,  in  a  mixed  integer  programming  problem 
or  a  nonlinear  programming  problem),  the  row/column  symmetry  of  the  Primal- 
Dual  .Algorithm  is  of  critical  importance  in  specializing  the  solution  approach.  The 
inherent  symmetry  of  the  algorithm  permits  easy  adaptation  to  branch-and-bound 
and  cutting-plane  approaches  to  mixed  integer  programming,  to  column  generation 
settings,  a.s  well  a';  to  primal  and  dual  decomposition  techniques. 
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We  believe  the  reason  for  this  flexibility  offered  by  the  algorithm  lies  in  its 
more  complete  mathematical  foundation.  There  is  a  natural  consistency  that 
arises  from  the  choice  of  a  vector  space  having  the  same  dimension  as  the 
problem  variables  that  is  lacking  in  other  approaches.  A  natural  geometric 
interpretation  of  the  solution  trajectory  follows  directly  from  this  development.  In¬ 
cidental  issues  such  as  finding  an  initial  basic  feasible  solution  and  dealing 
with  degeneracy  are  resolved  constructively  in  this  mathematical  frame¬ 
work.  Other  approaches  resort  to  unnecessarily  complicated  tangential  efforts. 

All  the  research  results  reported  here  can  be  developed,  with  some  effort,  in 
the  framework  of  the  classical  Simplex  Method.  However,  we  choose  to  present  these 
results  in  the  manner  of  their  development  -  the  mutual  Primal-Dual  view  presented 
by  Graves. 


i: 


III.  IMPLEMENTATION  DESIGN  OVERVIEW 


A.  INTRODUCTION 

We  seek  to  demonstrate  the  efficiency  of  row  factorization  (in  particular,  using 
dynamic-dimension  bases).  Accordingly,  we  will  implement  the  ideas  developed  here 
and  extensively  test  them  within  a  commercial- qu ah ly  optimization  system:  the  X- 
System  [1975]  employs  the  Graves  mutual  primal-dual  algorithm  in  a  variety  of 
large  scale  optimization  applications,  including  linear,  nonfinear,  mixed  integer  and 
decomposed  models.  The  benchmark  test  suite  is  drawn  from  a  wide  variety  of  actual 
applications,  and  our  goal  is  to  improve  the  efficiency  of  a,a  already  well-known  and 
highly  regarded  system. 

B.  DESIGN  CONSIDERATIONS 

W  e  want  to  test  our  ideas  by  repeatedly  solving  many  medium-  to  large-size  lin¬ 
ear  programming  problems  (i.e..  approximately  8.000-10,000  constraints  and  15.000- 
20.000  structural  variables).  Larger  problems  are  of  interest,  but  for  purposes  of  this 
research,  we  are  limited  to  a  relatively  modest  computer  budget  on  an  IBM  3033/.^? 
computer  under  the  \'M/CMS  operating  system  using  VS  Fortran  1.4.1.  We  wish 
to  support  the  computational  enhancements  common  in  commercial  mathematical 
programming  systems  (e.g..  bounded  variables,  ranged  constraints,  parametric  pro¬ 
gramming.  etc.).  \\'c  require  a  primal-dual  implementation  that  offers  complete 
flexibility  in  determining  solution  strategy.  In  addition,  each  experimental  imple¬ 
mentation  musl  support  all  the  routine  housekeeping  of  the  optimization  system 
(e.g.,  re-inversion,  crash  starts,  relaxations,  restrictions,  extrinsic  and  intrinsic  enu¬ 
meration.  etc.  '.  f-ina!l\.  we  seek  the  capability  to  identify  desired  ro\^•  factorization 
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structure  within  the  LP  model  instance,  either  by  communication  from  the  modeler 
or  automatically. 

C.  DESIGN  TEMPLATE 

To  establish  a  conceptual  framework  for  the  evaluation  of  our  algorithms,  it 
is  useful  to  outline  the  important  aspects  of  our  implementation  and  identify  the 
crucial  steps  which  most  strongly  influence  performance.  Recall  from  Chapter  2  our 
basic  tableau  in  partitioned  form; 


U)  Uj)  Ujj) 


(0 

A^i  Ai2 

(n) 
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where  we  have  made  the  substitution  —  B  */  in  the  right-hand  side  column: 
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The  data  requirements  of  the  algorithm  are: 


.3. 


Access  to  the  original  problem  data: 

.•\  representation  of  that  part  of  the  current  tableau  we  have  chosen  to  represent 
implicitly,  and: 

.A  represent  at  ion  of  t  hat  part  of  the  current  tableau  we  have  chosen  to  represent 
explicit  ly. 


We  new  (eti'ide:  f  ,-'].  ef  ;iie-e  rc'cpiireinent s  in  greater  detail, 
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Primal  simplex  implementations  typically  are  column-oriented  and  thus  re¬ 
quire  column- wise  access  to  the  problem  data.  In  the  primaJ-dual  method,  restric¬ 
tion  to  either  column-wise  or  row-wise  access  alone  exacts  a  serious  computational 
penalty.  Thus,  our  design  allows  both  column-wise  and  row-wise  access. 

The  design  of  an  efficient  large-scale  tableau  element  generation  representation 
is  remarkably  complex.  Many  subtle  software  engineering  and  hardware  environment 
issues  can  have  a  profound  influence  on  the  intricacy  2md  performance  of  promising 
designs.  The  design  excursions  reported  here  are  inexorably  influenced  by  architec¬ 
ture  of  the  host  computer,  operating  system  and  implementation  language.  However, 
the  reported  design  philosophy  has  been  tempered  with  experience  on  many  other 
computers  of  widely  varied  designs.  The  proposed  innovations  adapt  quite  well  to 
floating-point  pipelines,  large  cache  memories  and  parallel  architectures. 

Because  of  its  fundamental  role  in  the  construction  of  the  tableau,  we  maintain 
an  explicit  representation  of  /IJ"/  (we  thus  refer  to  v4i)  as  the  “explicit  kernel"). 
Although  any  of  the  various  techniques  such  as  LLKLDL^  or  QR  decomposition 
or  product  form  inverse  (e.g..  Golub  and  Van  Loan  [1985]  or  Magnanti  [197G]).  are 
suitable  for  representing  An  or  ,  a  difficulty  arises  in  this  algorithmic  setting. 
While  in  the  primal  Simplex  method  all  updates  to  the  basis  take  the  form  of  rank 
one  column  exchanges,  our  sotting  admits  more  general  updates,  which  include  single 
row  exchanges,  single  row  and  column  exchanges,  single  row  and  column  deletions, 
and  multiple  row  exchanges  of  Aj"]’  (multiple  columr  exchanges  of  An).  While 
this  does  not  preclude  the  use  of  any  particular  representation,  it  adds  a  level  of 
complexity  not  usually  encountered  in  more  traditional  implementations. 

We  also  maintain  an  explicit  representation  of  the  right-hand  side  column  and 
the  bottom  (cost )  row.  Because  of  the  symmetric  nature  of  the  mutual  primal-dual 
method,  a  sensible  apjiroaeh  is  to  al'cjcate  a  sinch'  storage  array  for  both  quantities. 


which  taken  together  are  called  the  “problem  rim”.  If  we  generate  an  explicit  rep¬ 
resentation  of  the  pi\-ot  row  and  pivot  column  of  the  tableau  at  each  iteration,  then 
we  may  update  the  problem  rim  array  using  the  simple  pivot  transformation.  By 
adopting  the  convention  of  labeling  the  nonbasic  constraints  as  rows  1  through  m 
and  labeling  the  basic  constraints  as  rows  (m  -|-  1)  through  (m  -|-  n)  (assuming  our 
problem  instance  has  m  structural  constraints  and  n  nonnegativity  constraints),  the 
problem  rim  array  is  partitioned  as  follows: 

1.  The  first  portion  of  the  array  holds  values  corresponding  to  the  nonbasic  con¬ 
straints  in  region  (z)  of  the  tableau  (3.1).  Since  these  are  nonnegativity  con¬ 
straints  and  they  are  nonbinding  (nonbasic),  the  values  in  region  (z)  of  the  rim 
array  are  those  of  the  currently  (possibly)  nonzero  variables. 

2.  The  second  portion  of  the  array  holds  values  corresponding  to  the  nonbasic 
structural  constraints  in  region  (it),  and  thus  the  values  in  the  rim  array  are 
the  current  slack  or  violation  in  these  constraints. 

3.  The  third  region  of  the  rim  arra%,  beginning  at  position  m  +  1.  corresponds 
to  basic  (binding)  structural  constraints,  and  thus  the  values  are  those  of  the 
corresponding  (po.^5ibly)  nonzero  dual  variables. 

4.  The  fourth  region  of  the  rim  array  corresponds  to  basic  nonnegativity  con¬ 
straints.  and  thus  the  values  arc  those  of  the  corresponding  (possibly)  nonzero 
dual  variables,  conventionally  called  “reduced  costs". 

The  re't  of  the  taleleau  we  represent  implicitly,  by  simply  recording  in  a  stor¬ 
age  array  the  currem  ordering  of  nonbasic  and  basic  constraints.  When  additional 


information  from  the  tableau  (for  example,  a  row  or  column)  is  required,  we  con¬ 
struct  it  from  our  representation  of  the  current  row  ordering  and  the  original 
problem  data. 

An  overview  of  the  solution  process  is  as  follows: 

1.  Identify  an  initial  basic  solution.  As  stated  previously,  the  origin  is  always 
such  a  solution,  and  thus  we  may  always  begin  with: 

5°  =  -I 
=  A, 

or  any  other  suitable  basis. 

2.  Check  for  optimality.  If  there  are  currently  no  primal  violations  {  a  negative 
value  in  the  first  or  second  region  of  the  rim  array)  and  there  are  no  dual 
violations  (a  positive  value  in  the  third  or  fourth  region  of  the  rim  array), 
then  the  current  solution  is  optimal.  Otherwise,  we  proceed  to  step  3. 

3.  Select  either  a  primal  or  a  dual  violation,  perform  a  pivot  which  makes  progress 
towards  reduring  that  violation  and  return  to  step  2. 

Since  we  desire  to  maintain  current  information  in  the  rim  array  by  means  of 
the  pivot  transformation  updates,  we  require  a  representation  of  the  pivot  row 
and  pivot  column  at  each  iteration.  Since  we  explicitly  maintain  only  Aj"/  and 
the  problem  rim.  we  see  that  a  key  computational  step  in  our  implementation 
is  the  generation  of  tableau  rows  and  columns. 

Recall  from  C'iiapter  2  the  principal  part  of  the  primal  tableau  (a  symmetric 
de\-elopmeiit  froii:  thf  dual  perspective  is  also  possible,  and  is  of  course  equi\'alent ): 
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where  P  =  —B~^  is  the  conjugate  row  basis, 
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and  D  contains  the  nonbasic  rows 


D=  /.'!  J  ^  (3.4) 

(jl)  >421  A22  J  • 

D.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 
EXPRESSIONS  IN  COLUMN  GENERATION 

Now  consider  the  generation  of  column  $  from  (3.2).  Rewriting  (3.2)  in  a 
manner  that  highlights  our  intentions; 
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By  properly  sequencing  our  computations  we  will  exploit  the  fact  that  region 

{ii)  of  a  given  column  is  simply  a  linear  combination  of  region  (?)  of  the  same  column. 

Assume  we  want  to  place  the  current  representation  of  column  s  into  a  work 

arra}'  r.  which  we  partition  as  =  (^1.22)^  to  correspond  to  (3.5).  Expressed  in 

terms  of  an  explicit  transformation  kernel  .47/.  we  first  compute  region  (?)  of  column 

s  as; 


if  s  is  ill  {]). 


or 


zi  =  if  5  is  in  (jj). 

Having  done  this,  we  then  compute  22  as: 


22  =  -A21Z1 


if  s  is  in  (j). 


or 


22  =  -A21Z1  +  (A22y  if  5  is  in  (jj). 

Assuming  an  LU  representation  of  An,  we  first  compute  region  (i)  of  column 


s  as: 


or 


LUzi  =  e’  if  column  s  is  in  (;), 


Li'zi  =  (An)''  if  column  s  is  in  (jj), 

Having  done  this,  we  then  compute  22  as: 


or 


Z2  —  — A2121 


if  column  s  is  in  {j), 


Z2  =  -A2121  +  (A22)’  if  column  s  is  in  [jj). 

Then  the  current  representation  of  column  5  is  available  in  z^  —  (21.22)^. 

■IS 


E.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 
EXPRESSIONS  IN  ROW  GENERATION 

The  computation  of  row  t  of  the  tableau  proceeds  in  a  similar  manner.  We 
now  view  the  principal  part  of  the  tableau  as: 


0) 


Ui) 


DP  = 


(0 

(ii) 


/ 


All  (^11^  )Ai2 


i-n 

^11 

1-1' 


\ 


If  we  want  to  place  the  current  representation  of  row  <  in  a  work  array  z 
partitioned  as  2  =  (23,24),  we  may  first  compute  region  (j)  of  row  t  as: 


53  =  (^4-’), 


if  row  i  is  in  (?), 


or 


23  =  (  — .42i)M7i*  if  row  t  is  in  (ii). 

We  then  compute: 


^3 


(-4,2) 


if  row  i  is  in  (i), 


or 


£4  i3(.4i2) -t- (.422)(  if  ro"’  ^  is  in  (ii). 

•■Mternately  u«inc  an  U'  representation  of  .4]]. 
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zzLU  =  tt 


if  row  t  is  in  (t), 


or 


zzLU  =  (— j42i)t  if  row  t  is  in  (ii). 

We  then  compute  24  as: 


^4  —  ■23(.4i2) 


if  row  i  is  in  (?), 


or 


24  =  53('4i2) -f  if  row  t  is  in  (iz), 

and  the  current  representation  of  row  t  is  available  in  2  =  (23,24). 

We  see  that  in  each  case  calculations  proceed  by  first  using  a  representation  of 
All  to  compute  a  portion  of  the  row  or  column  and  then  using  this  initial  computa¬ 
tion  and  original  problem  data  to  compute  the  remaining  part.  We  will  discover  that 
our  specializations  extend  this  approach  by  introducing  additional  tableau  partitions 
which  allow  this  computational  strategy  to  be  applied  on  a  larger  scale. 

As  previously  mentioned,  an  important  implementation  challenge  in  this  algo¬ 
rithmic  setting  is  the  dynamic  behavior  of  Au.  We  see  from  the  primal  row  basis 
(3.3)  and  nonbasic  rows  (3.4)  that  the  dimension  of  An  corresponds  to  the  number 
of  basic  structural  constraints,  or.  equivalently,  to  the  number  of  nonbasic  nonneg¬ 
ativity  constraint. (recall  that  if  a  nonnegativity  constraint  is  nonbasic  and  thus 


.bO 


nonbinding,  the  corresponding  variable  may  possibly  be  nonzero).  Recalling  that 
our  primal  view  of  a  pivot  is  as  an  exchange  of  constraints  between  B  and  D,  we 
see  that  one  of  tour  cases  may  occur  during  a  pivot: 

1.  A  structural  constraint  enters  the  basis  B  and  a  structural  constraint  leaves 
the  basis  and  enters  D.  Since  the  number  of  basic  structural  constraints  (and 
the  number  of  nonbasic  nonnegativity  constraints)  rentjains  unchanged,  the 
dimension  of  An  is  unchanged.  A  pivot  of  this  type  involves  a  row’  in  region 
(j)  of  B  (3.3)  and  a  row  in  region  (ii)  of  D  (3.4),  and  thus  it  corresponds  to 
a  pivot  coordinate  in  the  location  ((n),  (j))  of  the  tableau  (3.5). 

2.  A  nonnegativity  constraint  enters  the  basis  and  a  nonnegativity  constraint 
leaves  the  basis.  Again,  the  dimension  of  An  remains  unchanged.  Since  this 
pivot  involves  a  row  in  region  (jj)  of  (3.3)  and  a  row  in  region  (i)  of  (3.4),  the 
corresponding  tableau  (3.4)  pivot  coordinate  lies  in  ((0>(ii))’ 

3.  A  structural  constraint  enters  the  basis  and  a  nonnegativity  constraint  leaves 
the  basis,  and  thus  the  number  of  basic  structural  constraints  (equivalently, 
the  number  of  nonbasic  nonnegativity  constraints)  increases  by  one.  The  di¬ 
mension  of  .4]]  is  increased  by  one.  This  corresponds  to  a  pivot  coordinate  in 
region  {(ii).(jj))  of  the  tableau  (3.5). 

•1.  .4  iionnegativity  constraint  enters  the  basis  and  a  structural  constraint  leaves 
the  basis,  and  thus  the  dimension  of  An  is  decreased  by  one.  The  correspond¬ 
ing  pivot  coordinate’  in  (3.5)  is  ((0-(i))- 

W  e  see  tliat  we  may  exert  some  influence  on  the  behavior  of  the  dimension  of 
All  by  our  strategy  for  selecting  target  violations  for  primal  and  dual  constraints 
(i.e..  onr  j'riciiic  strategy!  and  through  our  tie-breaking  rules  for  choosing  pivot 


row/columns,  and  that  this  dynamism  is  an  inherent  feature  of  our  algorithm.  We 
have  already  seen  the  fundamental  importance  of  the  kernel  (/lu)  in  our  compu¬ 
tations.  Thus,  a  successful  implementation  must  manage  this  dynamic  behavior 
efficiently  and  reliably. 

To  illustrate  in  a  familiar  setting  the  challenge  this  offers,  consider  a  LU  rep¬ 
resentation  of  the  kernel  Au' 

=  LU. 

A  pivot  in  tableau  coordinate  ((?),  (jj))  results  in  a  column  exchange  in  An 
(3.6)  in  the  more  convenient  form; 

L-^An=U.  (3.7) 

When  column  re[)lace5  the  colutnn  of  /In  to  form  /In-  "'o  have 


where  q  =  has  replaced  the  column  of  I' .  We  now  must  restore  the 

ui)j>er  triangular  funn  of  I' .  We  would  prefer  a  method  which  demonstrates  strong 
nuriifTical  .stahili’’.  afid  v.hi'h  also  preserves  the  sparsity  of  I  .  S('\eral  methods 


(3.6) 
.  Writing 


have  been  proposed  (Bartels  and  Golub  [1969],  Forrest  and  Tomli".  [1972],  Saunders 
[1976]  and  Reid  [1982]).  Perhaps  the  most  widely  used  method  is  that  of  Saunders. 

Assume  for  the  moment  that  An  has  a  block  (or  “bump”)  triangular  form  (we 
discuss  how  this  is  done  shortly).  Then  An  appears  as  shown  in  Figure  3.1. 


Figure  3.1:  Bump  Triangular  Form  of  An 

Eacli  bump  consists  of  (possibly)  several  columns  that  extend  above  the  main 
diagonal.  Tliese  columns  are  called  “spikes'  .  and  the  tallest  spike  within  a  bump  is 
placed  in  the  right-most  column,  so  that  a  bump  appears  as  shown  in  Figure  3.2. 

The  block  triangular  form  of  An  results  in  a  LU  decomposition  which  has  the 
form  shown  in  Figure  3.3. 

Saunders  exploits  the  form  of  these  LV  factors  by  maintaining  a  permutation 
of  the  columns  and  rows  of  ('  so  that  all  the  spikes  appear  on  the  right-hand  side 


Figure  3.2:  Spikes  within  a  Bump 


0  =  p'^rp  = 


I  R  1 
0  F 


R 


F 


When  tlie  p'^‘  column  of  Au  replaced  by  ,  the  method  proceeds  as  follows; 

1.  Delete  the  p"  column  of  xF.  and  move  all  the  following  columns  one 
position  to  the  h  ft . 


2.  Place  o  =  /  ~'i’A  in  the  riclit -hand  column  of  I 


in 


A, 


■^1 


I  -i-  1, 


u 


Figure  3.3:  LV  Decomposition  of  Ax\ 

3.  Move  the  /)'*'  row  of  l\  ity.  to  the  ijottoin  of  A',  and  move  all  rows  in  between 
up  one  position.  Note  that  this  jneserves  the  triangular  form  of  all  but  the 
last  column  of  A,  and  ftirtlicr  note  that  row  Up  lias  nonzeros  only  in  columns 
corresponding  to  /A 

4.  Using  Gaussian  elimination  with  row  interchanges,  transform  i/p  to  a  sin¬ 
gleion  column,  thereby  restoring  the  u})[>et  triangular  form  of  I’.  Thus, 
f '  —  Er  ■  .  .  E]l  \  where  I  is  the  final  updated  form  of  E  and  £V.  •  •  .  Fi  file’  the 
elementary  transformation  matrices  corresponding  to  the  Gaussian  elimination 
steps  (s(’<>,  e.g.,  Murtagh  [lf>8lj). 

5.  -Apply  these  saim’  (iaussian  elimination  steps  to  A”',  forming 


/  ‘  - 


/....A,/.-'. 


Typically,  JL"*  is  held  in  product  form  and  thus  these  transformations  are 
simply  added  to  the  current  list  of  transformation  vectors  (“eta”  vectors). 

This  method  hais  the  virtue  that  new  nonzeros  can  be  created  only  in  the 
submatrix  F  of  U ^  and  thus  R  may  be  carried  in  peripheral  storage.  The  numerical 
properties  are  reported  to  be  quite  good. 

The  computational  burden  of  continuously  maintaining  A\\  in  block  triangular 
form  is  excessive,  and  the  usual  approach  is  to  refresh  the  structure  as  part  of  the 
bzisis  re- inversion  routine.  Between  the  re- inversions,  the  structure  is  left  untended. 
An  effective  heuristic  due  to  Hellerman  and  Rarick  [1972]  is  commonly  used  for  this 
purpose. 

A  pivot  in  tableau  coordinate  ((n),(i))  results  in  a  row  exchange  in  A\\. 
Writing  (3.6)  as 

AuF-'=L  ■  (3.8) 

VMien  low  replaces  the  row  of  An  to  form  ,4]].  have 


i  L 

I  _ 

where  row  A  =  replarrs  the  q’^‘  row  of  /,.  The  desired  structure  of  L  may  be 

restored  by  meiln'ds  syini!)elric  to  tho.'^e  cicveloj>ed  for  the  column  e.xchange  case, 
again  formiim  : 


An  =  10  . 

A  pivot  in  tableau  coordinate  causes  the  dimension  of  An  to  in¬ 

crease  by  one  through  the  addition  of  a  row  and  a  column.  It  is  convenient  to  add 
the  new  row  to  the  bottom  of  An  and  the  new  column  to  its  right.  If  An  is  the 
new  kernel,  then 


An  — 


An  fl* 


where  0^,0*'  and  ar,k  are  original  matrix  coefficients. 

The  desired  updated  decomposition  is  of  the  form; 


(3.9) 


(a  ' 

An  a 

'  L  0  ' 

■  U  u* 

^r,k 

Ir  Ir.r  , 

0  Uk,k  _ 

which  requires  that 


LiA  = 


IrO  =  Or. 


and 


Or^k  —  "|-  Ij- 


■Setting  =  1  -  we  haxi  : 


The  final  pivot  coordinate,  ((?)i(j))  results  in  the  dimension  of  An  decreasing 
by  one.  If  An  is  the  current  transformation  kernel  and  An  results  after  the  pivot, 


without  loss  of  generality  we  have: 


All  — 


0.r,k  0-T 

h,k  0 

Ur.r  Ur 

irr 

I - 

1 _ 

0  Ui 

(3.10) 


where  (Or.it,  Or)  is  the  leaving  row,  (or.it,  the  leaving  column  and  ar,k  the 

pivot  element.  Using  an  analysis  similar  to  Rosen  [19C0],  we  note  that  if  the  square 
matrix  C  is  nonsingular  and  is  partitioned  as: 


^  _  f  c,  C2  ■ 

[  C3  C,\  , 

where  Ci  is  square  and  nonsingular.  C4  is  square  and  Cq  —  C4  —  C3Cf'C2  is  non¬ 
singular.  then 

_  ■  cf’  +  cr'C2Co-'C3cr’  -cr'GCo ' ' 

“  -Co’CsCr’  Co’  _  . 

.Applying  this  re.^ult  to  (3.10).  we  discover  that 

--111  -  O  (Ut-Jc)  Or  “  Z/  1  f  1 


or 


If  tliV 
hand . 


terin 

Unfortuiia'c 
Ir.'  jx'rfi jr'i i;n 


cr, 


k 

G  Ur 


(3.11 ) 


should  be  tlie  zero  matrix,  our  new  decomposition  is  at 
need  not  be  zero,  but  v.c  may  guarantee  it  to 
r  a  jiT' ii:.'i!:.ar\'  rolunui  ('\c}iaiie('  update  tri  .-I,;.  We  rei>la' - 


column  ^ar,t,a*)  with  column  (1,0)^  and  compute  the  updated  factors  exactly  as 
we  did  in  the  column  exchange  case  considered  earlier.  This  results  in  a  modified 
transformation  kernel  An  with  factors 

2  _  1  flr  _  'lk,k  0  Ur,r 

"  "  [  0  iu  J  “  [  4  J  [  0  Lh  . 

Then  the  second  term  of  (3.11)  is 

=0 

and  thus  the  updated  LU  factors  are  at  hand: 

An  =  L\l  \  ■ 

W’e  now  sec  the  advantages  and  disadvantages  of  the  Graves  mutual  primal- 
duai  method.  It  has  the  advantage  that,  although  the  transformation  kernel  can 
be  as  large  as  rn.  the  number  of  structural  constraints  in  the  problem  instance,  its 
dimension  is  actually  equal  to  the  current  number  of  binding  structural  constrain’ ' 
(or  equivalently,  to  the  number  of  potentially  nonzero  primal  variables).  .At  our 
initial  basic  solution,  the  kernel  dimension  is  zero.  If  the  maximum  kernel  dimen¬ 
sion  during  the  course  of  the  solution  trajectory  is  much  smaller  than  rn.  we  enjoy 
significatit  comjiut at iona!  advantages  in  storage.  ui)date  and  comjjutation.  For  a 
great  man\'  I.l’  mod*'!  instances,  this  is  precisely  the  case. 

The  dynamif  nature  of  the  transformation  kernel  clearly  complicates  update 
[)ro(  edures.  Whih'  th'  usual  primal  S.mplex  method  requires  only  th('  column  ex- 
rhaniT''  ujidaie.  (ira\f  -  iii'’thi>d  require-  fi.nr  update  case-. 


In  our  implementation,  we  will  seek  methods  for  handling  the  transformation 
kernel  that  enjoy  the  advantages  while  mitigating  the  disadvantages. 

It  is  this  general  strategy  for  row’  and  column  generation  we  develop  and  en¬ 
hance  in  our  specializations.  By  identifying  a  special  structure  in  the  problem  data, 
we  will  introduce  additional  hnear  dependencies  in  the  rows  and  columns  of  the 
tableau  which  will  further  simplify  their  generation  and  will  also  further  reduce  the 
dimension  of  the  transformation  kernel. 


IV.  FACTORIZATION 


A.  INTRODUCTION 

Each  of  the  algorithms  we  present  in  this  research  can  be  viewed  as  a  special¬ 
ization  of  a  general  approach  to  large-scale  linear  programming  developed  by  Graves 
and  McBride  [1976],  which  they  call  the  “factorization  approach”.  It  is  based  on 
distinguishing  special  rows  and  columns  in  a  way  that  allows  large  parts  of  the  ba¬ 
sic  tableau  to  be  represented  implicitly,  to  be  generated  easily  from  the  remaining 
explicit  parts  only  when  actually  required  by  the  algorithm. 

Although  sometimes  used  interchangeably  in  the  literature,  we  recognize  a 
distinction  between  partitioning  methods  and  factorization  methods.  In  the  former, 
a  formal  distinction  is  made  between  substructures  which  appear  in  the  model  in¬ 
stance,  usually  constraints  or  variables.  .An  approach  for  solving  the  problem  is  then 
developed  wliich  exploits  the  substructures.  For  example,  Dantzig-V’olfe  decompo¬ 
sition  is  such  an  approach  which  partitions  constraints,  while  in  Benders  decompo¬ 
sition  the  variables  are  partitioned.  Such  an  approach  may  be  applied  statically,  in 
which  case  the  desired  partitions  are  identified  at  the  start  of  the  solution  process 
and  remain  fixed  throughout,  or  it  may  be  dynamic,  in  which  ca'^e  the  partitioning 
may  be  adjusted  as  the  solution  proceeds. 

In  contrast,  we  consider  factorization  methods  to  be  those  in  which  a  formal 
distinction  is  made  between  substructures  in  the  LP  ba^is  (and  thus  in  the  basic 
tableau),  The  aleorithrn  is  then  specialized  to  exploit  this  special  substructure. 
1  hu:-.  w.  cie\>  tie'  Gl'P  algorithm  of  Daiiizig  and  \’an  .Slyke  [1967]  and  the  mul- 
1 1( ommodii y  netweik  algisriihm  of  Hartman  and  Lasdon  [1970:  as  example-s 


of  factorization  algorithms.  Factorization  methods  may  also  be  static  or  dynamic, 
depending  on  whether  the  dimension  of  the  factored  substructures  in  the  LP  basis 
may  vary  during  the  solution  process. 

Our  research  is  concerned  with  dynamic  row  factorization  algorithms.  We 
develop  the  general  dynamic  factorization  approach  here  and  discuss  the  important 
aspects  of  the  algorithm.  In  subsequent  chapters,  we  specialize  this  general  approach 
to  each  of  several  special  LP  row  structures.  This  general  development  closely 
parallels  that  of  Graves  and  McBride  [1976]. 

B.  THE  FACTORED  TABLEAU 

The  problem  to  be  considered  is: 

(FLP)  min  :  wy 

s.t.  :  Ey  <  r  }  explicit  constraints 

Ey  <  b  }  factored  constraints 

<  0  )  nonnegativity  constraints  , 

where  y  is  a  ??  vector  of  decision  variables,  u-  is  a  n  vector  of  cost  coefficients.  E 
is  an  m  by  n  matrix  of  constraint  coefficients  for  ‘"explicit"  constraints  with  right- 
hand  side  m  vector  r.  F  is  a  p  by  n  matrix  of  constraint  coefficients  for  “factored" 
constraints  with  right-hand  side  p  vector  b.  and  —I  is  the  negative  of  the  n  by  n 
identity  matrix.  In  this  general  development,  we  refer  to  the  F-type  constraints 
as  “factored"  oiil\'  to  distinguish  them  from  the  ‘"explicit"  F-tvpe  constraints,  and 
assume  nothing  about  their  structure.  Not  until  our  specializations  in  later  chapters 
will  we  imf)Ose  special  structure  on  F,  and  the  structures  we  will  consider  permit 
the  represeiUat ion  of  the  F  type  constraints  without  the  inversion  of  a  matrix.  We 
will  see  that  this  apfiroach  is  centered  around  handling  the  part  of  the  basis  cor¬ 
responding  to  the  F-type  constraints  exjdicitly  while  fartorinc  the  portioii  of  tlie 

(ij 


b2Lsis  corresponding  to  the  F-type  constraints.  The  notation  is  chosen  to  suggest 
this  idea. 


As  we  saw  in  Chapter  2,  in  the  mutual  primal-dual  method  the  primal  al¬ 
gorithm  and  the  dual  algorithm  share  the  same  tableau,  and  thus  the  tableau  for 
(FLP)  may  be  derived  from  either  perspective.  We  present  only  the  derivation  from 
the  primal  viewpoint  here. 

Recall  that  a  basis  for  the  primal  algorithm  consists  of  n  linearly  independent 
row's  from  the  constraint  coefficient  matrix  when  it  is  assumed  to  include  both 
structural  (explicit  and  factored)  and  nonnegativity  constraints.  Assume  that  the 
current  row  basis  consists  of  /  rows  from  F,  k  row's  from  F  and  {n  —  k  —  1)  rows 
from  —I.  Using  the  notation  of  Chapter  2  temporarily,  we  have; 


i+k 

n~l—k 

( 

\ 

B  = 

An 

A]2 

}  1  k 

\  0 

J 

}Tt  -  1  —  k 

where  [,4ii  /I12]  includes  al!  basic  structural  rows,  from  both  E  and  F. 

We  will  ultimately  be  interested  in  isolating  the  effect  of  each  type  of  struc¬ 
tural  constraint  algebraically  in  the  factored  tableau,  and  thus  we  require  greater 
resolution  in  our  factored  basis.  Introducing  obvious  notation,  we  have: 


wh^'ie 


i.du  .-1 


1-.. 


k 
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Eu 


r.-t~k 
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E\-\  / 
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:.4 


)  1 , 


f  - 11  f- 1 

f  1 1  } 
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From  the  development  in  Chapter  2,  we  recognize  An  as  the  exphcit  kernel, 
and  thus  exists.  It  then  follows  that  it  is  always  possible  to  identify  among  the 
columns  of  [Fn  F^]  a  nonsingular  submatrix  Fn  of  dimension  k,  since  otherwise 
the  rank  of  [Fn  F12]  is  at  most  {k  —  \  )  and  thus  the  rank  of  B  is  at  most  (n  —  1). 
We  will  later  see  that  one  of  the  important  implementation  challenges  is  the  task  of 
efficiently  managing  the  structure  and  nonsingularity  of  Fn. 

The  full  factored  row  basis  is  then: 


n—/— A: 


B  = 


Introducing  the  notation: 


U) 

En 

E\2 

E\z 

M 

j 

Uj) 

fn 

F,2 

Fn 

]k 

Ujj) 

[  0 

0 

-J  ) 

}  n  —  1  — 

(4.1) 


.4ii  —  Fi2  —  Fi]Fjj^Fi2 

.4,2  =  F,3-F„F-^F,3  , 

its  inverse  is: 

~-C]i'Fi2-4i/  (/ +  Fii’Fi2/4i/F]i)Fj,’  I' —  Fi2A-i^  An) 

_  A^\  -^1,'FnFf/  A^i  A\2 

0  0  -1 

I  Grouping  the  cocfiicients  of  the  nonbasic  constraints  and  applying  the  same 


column  ordering  s'iehJ': 


1-i-k 


D  = 
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-1 

0 
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F2X 

F22 

F23 

}p-k 

(iit) 

,  0 

-I 
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(tv) 
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E22 

E23  } 

]m  —  1 

(4.2) 


We  will  explain  the  ordering  of  the  rows  of  D  shortly. 

The  principal  part  of  the  factored  tableau  is  DP,  where  P  =  is  the 

conjugate  row  basis.  Introducing  the  additional  notation: 


F22 

A 

F22 

F23 

A 

^23 

-F2,F-^Fr^ 

A2] 

A 

E22 

-E2iFu'Fn 

A22 

A 

E23 

~  E2\F IX  F\3  , 

the  principal  part  of  the  factored  tableau  is: 

(.)>  (n) 


(;;j) 


(0 


-Fu'FuAu  {I  ^  F::FnA-,lEu)Fu^  F^,\F,,  -  FnA'^An) 
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(t?) 

~-^22-4ii‘ 

t  ^'22-4ii''^^ 

u-F2x)Fx-x^ 

P23  —  F22AXX  Ax2 

(?n) 

-■An' 

--■An 

■TnT.V 

^11' -^12 

(n') 

\ 

--'A21-'Aii' 

(MxA^^E 

11  —  E2X  )/■  11' 

A22  —  '421  Axx  Ax2 

j 

DP  = 


(4.3) 

Part  it  ioniiic  ic  =  ( ti'i .  le^.  ).  =  (r'l.rn)^  and  —  (^1-^2)^ 

int  rodiiring  t  he  not  at  iwn: 


li" 


»'■.  -  U  i  /  ,1  /  1. 


U'3 

A 

U.’3  - 

■  ty2F,,’F]3 

A 

b,~ 

FnF-^b, 

f'l 

A 

fl  - 

F„F,-,’6, 

^2 

A 

r2  ~ 

F2iF,^’6,  , 

the  complete  factored  tableau  is: 


(/  +  F,VF„.4-'£n)Fl' 

F 1i*(£i3  —  £i224],’,4i2) 

Fn](^  ~ 

~F22-^U 

( A2.4  r;  £„-£„)£, V 

£23  —  £22.i4j,'  /ti2 

-•-InF.F.-;' 

T)i'  y4i2 

^r.'A 

(.-l2..4r,'£,:  -£:,)£.V 

T22  ~ 

-U'2.4F 

(uj.-t),’  £]i  —  wi)F 11' 

W3  ~ 

(4.4) 


We  see  from  (  l.-l  l  that  witii  knowledge  of  the  current  factorization,  we  can 
construct  the  entire  tableau  from  F,”,’,  /If,’  and  the  original  problem  data.  The 
dimension  of  F,’,'  i.^-  eijual  to  the  number  of  F-type  constraints  that  are  currently 
basic,  and  thus  can  t,>e  at  most  p.  The  dimension  of  Tf,'  is  equal  to  the  number  of 
F-type  constraint''  that  are  currently  basic.  Hence,  its  dimension  cannot  exceed  m. 
We  call  ,4f,'  the  explicit  transformation  kernel  and  F,”,’  the  factored  transformation 
kernel. 

We  have  sf-en  in  ('hai)ter  that  the  origin  is  always  a  convenient  initial  basic 
solution  for  the  iiiulua!  i>rimal-dual  meihod.  and  the  same  is  true  for  the  factoriza¬ 
tion  aj){)roac!i.  An  imAa!  ba^ir  sohilitui  is  always: 


f.n 


=  [-/] 


Starting  from  this  solution,  we  may  view  the  algorithm  as  progressing  by  exchanging 
rows  of  F  and  E  from  D  with  nonnegativity  constraints  from  B.  We  will  find  it  useful 
to  associate  with  each  F-type  constraint  in  B  a  unique  nonnegativity  constraint  in 
D.  Borrowing  the  terminology  of  Dantzig  and  Van  Slyke  [1967],  we  will  call  each  such 
unique  noimegativity  constraint  (and  its  corresponding  variable)  the  “key”  variable 
for  the  associated  F-type  constraint.  The  algebraic  structure  arising  from  this  basic 
F-lype  constraint/nonbasic  “key”  variable  association  allows  us  to  represent  parts 
of  the  tableau  implicitly  rather  that  explicitly.  To  distinguish  these  “key”  variables 
in  D  from  others,  we  place  them  in  region  (a)  of  D.  This  is  the  reason  for  the  row 
ordering  mentioned  earlier. 

C.  BENEFITS  OF  FACTORIZATION 

Graves  and  McBride  [1976]  report  three  principal  benefits  of  the  factorization 
approach; 

1.  Good  starting  bases; 

2,  Reduction  in  memory  requirements; 

■3.  Reduction  in  work  per  pivot. 

.Althougli  the  origin  i.s  always  available  as  an  initial  basic  solution.  Grave.*;  and 
McRnde  [1976;  su  gcc'-t  that  a  better  one  can  frequently  be  found  with  only  small 


computational  effort  using  this  approach.  They  form  the  Lagrangian  Dual  of  (FLP) 
with  respect  to  the  E-type  constraints,  we  obtain: 

(LFLP)  mm  :  wy  -}-  A(Ej/  —  r) 

s.t.  :  Fy  <  h 

-/y  <0 

where  A  >  0  is  a  m  vector  of  dual  variables.  If  A  equals  A*,  the  optimal  dual  variables 
for  (FLP),  then  an  optimal  solution  of  (FLP)  is  among  the  optimal  solutions  of 
(LFLP).  They  speculate  that  if  A  is  a  “good”  estimate  of  A*,  then  using  A  in  (LFLP) 
and  solving  for  y  should  yield  a  “good”  starting  point  y  for  (FLP).  If  no  such  estimate 
A  is  available,  then  A  =  0  may  be  used.  Depending  upon  the  structure  of  the  E-type 
constraints,  the  resulting  problem  may  either  be  much  easier  to  solve  than  (FLP)  and 
thus  may  be  solved  optimally,  or  heuristics  may  be  used  to  find  a  computationally 
inexpensive  “good”  solution.  All  factorization  algorithms  can  be  started  by  solving 
(LFLP)  first. 

Standard  simplex  techniques  require  a  basis  of  dimension  (m  +  p)  to  solve 
(FLP).  The  factorization  approach  requires  the  explicit  transformation  kernel 
whose  dimension  is  at  most  m,  and  the  factored  kernel  Fu-  wnoso  dimension  is 
at  most  /).  Thus,  it  is  clear  that  we  expect  a  considerable  reduction  in  memory 
requirements  using  the  factorization  approach.  When  we  specialize  the  approach 
for  models  in  which  the  E-type  constraints  have  special  structure  (c.g..  GUB.  pure 
network  or  generalized  network),  we  will  see  that  memory  requirements  can  be 
further  reduced. 

W  hile  it  is  difficult  to  measure  the  computational  effort  per  pivot  when  product 
form  in\erse  or  basis  decomposition  techniques  are  used,  there  arc  indications  that 
the  factorization  ap[)roach  can  yield  substantial  benefits.  Both  analytic  (.McBride 


[1972])  and  empirical  (Tomlin  [1972])  studies  indicate  that  when  the  f-type  con¬ 
straints  are  GUB  constraints  (we  will  cover  this  in  detail  in  Chapter  6),  the  factoriza¬ 
tion  approach  results  in  per-pivot  computational  benefits  when  the  GUB  constraints 
comprise  at  least  30%  of  the  total  number  of  structural  constraints.  With  some  sim¬ 
plifying  assumptions,  a  similar  analysis  can  be  developed  for  the  case  where  the 
T-type  constraints  are  pure  network  rows. 

We  will  consider  each  of  three  different  factored  row  structures  in  this  work. 
The  complexity  of  the  factored  row  structure  determines  the  computational  difficulty 
of  representing  and  updating  the  factored  kernel  Fn.  The  simplest  row  structure 
we  consider,  in  which  the  factored  rows  are  generalized  upper  bounds,  allow’s  to 
be  interpreted  as  a  simple  permutation  matrix.  The  second  factorization,  in  which 
the  F-type  constraints  are  pure  network  rows,  is  a  more  complex  structure  which 
subsumes  generalized  upper  bounds  as  a  special  case.  In  this  case,  Fn  may  be 
interpreted  as  lepresenting  a  partial  ordering  defined  over  a  subset  of  the  rows  of  F. 
The  final  factorization,  in  which  the  rows  of  F  are  generalized  network  rows,  is  the 
most  general  we  consider  in  this  work  and  subsumes  both  the  other  factorizations. 
Here,  Fu  mav  be  interpreted  as  a  linear  transformation  in  which  a  near  partial 
ordering  exists  among  a  subset  of  the  rows  of  F. 


V.  GENERAL  IMPLEMENTATION  TOOLS 


A.  INTRODUCTION 

We  DOW  provide  an  overview  of  our  implementation  of  the  factorization  method. 
We  describe  the  representation  and  update  of  the  three  important  algorithm  compo¬ 
nents:  the  basic  tableau,  the  factored  kernel  and  the  explicit  transformation  kernel. 
The  fundamental  update  steps  are  described  in  terms  of  functions  which  operate 
'll  each  of  the  algorithm  components.  W'here  the  details  of  the  update  actions  are 
common  to  all  three  factorizations  treated  in  this  research,  the  details  are  presented 
in  this  chapter;  otherwise,  they  appear  in  subsequent  chapters. 

B.  THE  FACTORED  TABLEAU 

As  before,  .ve  are  inteicciod  in  tiic  prouiem. 

min  :  wv 

y 

s.t.  :  Ey  <  7-  }  explicit  constraints 

E  b  ]  factored  constraints 
~^y  E  0  }  nonnegativity  constraints 

Recall  the  algebraic  representation  of  an  arbitrary  primal  row  basis: 
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with  the  corresponding  nonbasic  rows: 
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F22 
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The  conjugate  row  basis  inverse  is: 


(5.2) 


(5.3) 


P  =  -B-' 
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and  tfic  j)riiicipal  pait  of  liic  fartoied  tablrau  is: 
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The  entire  principal  part  of  the  factored  tableau  can  be  constructed  with 
knowledge  of  the  row  and  column  ordering,  (or  Fu)  and  .  Our  discussion 
of  the  implementation  centers  around  the  representation  and  update  of  these  three 


components. 

We  restrict  our  attention  to  the  principal  part  of  the  tableau  because  our 
implementation  maintains  a  current  representation  of  the  right-hand  side  column 
and  bottom  (cost)  row  (collectively  referred  to  as  the  “problem  rim”)  at  all  times. 
These  quantities  are  always  immediately  available  and  need  not  be  computed  upon 
demand.  This  is  in  contrast  to  (5.4),  where  only  and  Fn  are  always  available 
and  any  other  quantity  of  intercut  (for  example,  a  complete  row  or  column)  must 
be  computed. 

C.  ALGORITHM  OVERVIEW 

In  broad  terms,  the  algorithm  proceeds  as  follows: 

1.  Establish  an  initial  basic  solution. 

2.  Stop  if  the  rurrciit  solution  is  optimal.  Optimalit}’  exists  when  the  right-hand 
side  colum.!  is  nonnegativ^’  and  the  bottom  (cost)  row  is  nonpositive. 

Based  oii  a  juicing  scheme  and  a  ratio  test,  select  and  generate  a  row  and 
colunin  of  i  .5. !  >. 

4.  Stop  if  thi'-  row  ann  column  rcvea!  infeasibility  or  unboundedness. 

o.  1  r -.n^f ’i::.  juwiu'-m  riu;.  tho  lab’eau  renresentation,  /n  and  .41]'.  and 
roi  uri,  ' ' '  s' '  ;  ■  J, 


'  i'.  ^  I 


ur'':-  ai.d  m'cessarv  Ujidate  operatio:,'.  we  will 


D.  IMPLEMENTATION  CONVENTIONS 


We  will  first  discuss  the  implementation  of  each  of  the  three  algorithm  com¬ 
ponents  individually,  and  then  provide  a  detailed  presentation  of  the  complete  al¬ 
gorithm.  Our  implementation  is  coded  in  the  FORTRAN  language,  and  thus  the 
following  discussion  is  colored  by  the  character  of  that  language. 

1.  The  Tableau 

The  mutual  primal-dual  method  fosters  a  symmetric  view  of  the  hnear 
programming  problem  and  the  solution  process.  This  symmetry  is  reflected  in  the 
indexing  conventions  required  in  our  implementation.  Assuming  now  that  (5.1) 
consists  of  m  structural  constraints  (both  factored  and  explicit)  and  n  structural 
variables,  we  require  the  structural  constraints  to  be  indexed  from  1  to  m,  and  the 
variables  from  (m  -f-  1)  to  {m  +  n).  While  our  statement  of  the  problem  specifies 
the  nonnegativity  constraints  explicitly,  we  will  of  course  deal  with  them  implicitly. 
Note,  however,  that  we  have  two  equivalent  interpretations  of  the  indices  (m  4- 
1)  through  (m  +  n).  \Vc  can  view  them  as  the  indices  of  the  structural  problem 
variables  or  as  the  indices  of  the  nonnegativity  constraints  corresponding  to  each 
of  the  structural  variables.  W'e  will  u.se  both  interpretations  in  our  development. 
Observe  also  that  wo  require  no  explicit  representation  of  logical  (slack,  surpdus  or 
artificial)  variables. 

The  indexing  of  the  original  problem  data  exists  independently  of  the 
solution  process  and  is  thus  referred  to  as  the  “extrinsic”  indexing.  The  factoring  of 
the  constraints  to  form  and  /)  re[)resents  a  second  indexing  of  tin'  [Problem  which 
defines  the  current  solution,  infiependent  of  the  extrinsic  indexing.  \\i.'  reh'r  to  thii- 
indexing  sy^lr’n,  a'  ‘'int  '  .  and  it  partially  defines  the  current  represent  at  lun  of 
If  and  I>.  Our  convention  fm  intrinsic  indexing  has  the  row^  of  />  labeled  from  1 
to  II:  and  1  i;e  ri j'.'.  -  of  }i  l<i : lei'-d  fi om  [ni  -  1  i  to  i  rn  ri  ■ . 


The  original  problem  data  is  stored  in  super-sparse  form,  (in  super-sparse 
form,  each  real  value  is  stored  once,  and  each  nonzero  coefficient  is  represented  by  a 
row  index,  a  column  index,  and  a  pointer  to  a  real  value)  with  unique  nonzero  real 
values  stored  once  in  an  array  of  type  DOUBLE  PRECISION,  and  referenced  by  the 
FORTRAN  equivalent  of  a  pointer  from  each  constraint  matrix  coefficient.  These 
pointers  to  the  real  values  are  singly-linked  by  both  row  (with  associated  column 
indices)  and  by  column  (with  associated  row  indices).  This  allows  symmetric  row¬ 
wise  and  column-wise  access  to  the  data  to  the  performed  efficiently,  which  is  an 
important  virtue  in  a  primal-dual  design. 

To  specify  the  structure  of  (5.4)  we  require  a  representation  of  the  cur¬ 
rent  mapping  between  the  extrinsic  indexing  of  the  problem  data  and  the  intrinsic 
indexing  of  the  current  tableau.  We  represent  this  mapping  by  two  arrays  of  type 
INTEGER,  each  dimensioned  from  1  to  (rn  +  n).  which  are  used  as  pointers.  We 
denote  these  arra\s  MI.N'T()  and  MEXT().  MI.NT()  represents  the  mapping  from 
intrinsic  to  extrinsic  indices,  and  MEXT()  represents  the  inverse  mapping.  Thus, 
if  ll.X'I  is  an  intrinsic  index  of  a  row  or  column  of  (5.4)  and  lEXT  is  the  extrinsic 
index  of  the  constraint  or  variable  currently  mapped  to  position  IINT  of  the  tableau, 
then: 

Ml.N’TdlNT,  =  lEXT 
MEXT(IEXT)  =  IINl 
Mi.\  ri.Mr.xT(ii:xr)i  =  lEXd 
MI  X'!  (MIN  Tdl.N  Tj!  =  IIN'I  . 

'1 ' '  f 'unpl'i''  till'  rcpri-'-ent  at  ion  of  (.5.1)  we  recpiire  knowledije  of  the  far- 
I  <  i:  ;,a :  i' ■:  riji:'"' nO'.;  alceliraieal;-.-  li\  repions  f/',  eti'.  We 

I  1 


handle  this  with  the  use  of  several  variables  of  type  INTEGER  which  are  used  as 
pointers  to  record  the  ending  intrinsic  index  of  the  various  tableau  regions.  Table 
(5.1)  lists  the  regions  and  the  associated  pointer  names. 

TABLE  5.1:  Indices  of  Tableau  Regions 
REGION  BEGINNING  INDEX  ENDING  INDEX 


(i) 

1 

MKC 

(ii) 

MKC  +  1 

MFR 

(iii) 

MFR  +  1 

MXR 

(iv) 

MXR  +  1 

M 

(j) 

M  +  1 

NXR 

(jj) 

NXR  +  1 

NFR 

(jjj) 

NFR  -f  1 

M  +  > 

The  arra\'.s  MI.\7'()  and  .MEXT{)  and  the  pointers  MKC.  MFR,  MXR. 
NXR  and  NFR  comprise  the  data  structures  necessary  to  represent  the  principal 
part  of  the  tableau.  7b  complete  the  representation  of  the  tableau  we  require  a  rep¬ 
resentation  of  the  right-hand  side  column,  the  bottom  row  and  the  current  value  of 
the  extremal  function.  The  right-hand  side  column  and  bottom  row  are  represented 
by  an  array  of  type  DOFBLE  PRECISION,  dimensioned  from  1  to  (m-fn).  referred 
to  as  RlMi  ).  Position^;  1  through  rn  correspond  to  the  right-hand  side  column  and 
posit  ion‘;  (rn  +  1  j  through  {iv  +  n; )  correspond  to  the  bottom  row.  Note  that  these 
are  intrinsic  inrlices.  7  he  curretit  value  of  the  extremal  function  is  held  in  a  scalar 
variable  of  type  DOl'BLE  PRECISION  called  OPT. 

W  e  now  discus-  the  operations  necessary  to  maintain  and  update  the 
ta'deau  represvntat  i>  i:.  7"  re  ognize  a'  z!  scdution  w('  simpl\'  scan  tlie  RlMi  j 

arra\.  \\  in'ii  th'-  ii.  p<i-itiuns  1  through  rn  are  nonnegat  i  \'i‘  and  tlmse  in 

prisi;  lull-  i  111  lit  hiiiuch  ( III  -r  7;  i  are  noiijKisit  im'.  t  he  current  solut  loti  i-  opt  imal. 


I  > 


Violations  of  either  sign  discipline  may  serve  as  potential  primal  algorithm  target 
rows  or  dual  algorithm  target  columns,  respectively. 

Having  selected  a  target  row  or  column,  we  must  generate  the  correspond¬ 
ing  row  or  column  of  the  tableau,  compute  ratios  with  respect  to  the  right-hand  side 
column  or  bottom  row’,  and  thus  identify  the  index  of  a  corresponding  column  or 
row.  We  then  generate  that  column  or  row  of  the  tableau.  Because  the  details  of  row 
and  column  generation  differ  according  to  the  factorization,  we  defer  their  discussion 
until  the  appropriate  chapter.  To  represent  the  update  operations,  however,  we  de¬ 
fine  the  abstract  functions  Generate_Row(IPRI)  and  Generate_CoIumn(IPCI). 
These  functions  accept  as  arguments  the  intrinsic  row  (IPRI)  or  column  (IPCI) 
index  of  the  current  tableau  and  return  the  current  representation  of  that  row  or 
column.  We  require  a  working  array  of  type  DOUBLE  PRECISION  to  hold  these 
current  representations.  We  define  ROWCOL()  to  be  such  an  array,  dimensioned 
from  1  to  (m  +  n).  where  positions  1  through  m  correspond  to  column  IPCI  of  (5.4) 
and  positions  (m  -r  1)  through  [m  -f  n)  correspond  to  row  IPRI.  ROWCOL()  is 
indexed  intrinsically  in  conformance  with  our  convention  for  (5.4). 

Recall  that  our  interpretation  of  a  pivot  from  the  primal  algorithm  view¬ 
point  is  as  an  exchange  of  rows  between  B  and  D.  The  row  indices  of  D  correspond 
to  the  row  indice.'  of  (.o.  l  i,  while  the  row  indices  of  B  correspond  to  column  indices 
of  (5.4  ).  Thus,  a  pivot  requiring  the  exchange  of  a  row  of  B  and  a  row  of  D  inquires 
the  exchange  of  a  row  index  of  (5.4)  with  a  column  index  of  (5.4).  Using  MINT() 
and  ME.\T().  such  an  exchange  requires  just  a  few  a.ssignment  statements.  For 
example,  if  the  pivot  rcr.v  IPHI  corresponds  to  extrinsic  index  EPRl  and  the  pivot 
column  IPCI  corrc'poini'  to  extrinsic  index  EPCl.  then  prior  to  a  pivot  we  have: 


MINT(IPRI)  =  EPRI 
MEXT(EPRI)  =IPRI 
MINT(IPCI)  =  EPCI 
MEXT(EPCI)  =IPCI 


and  after  the  pivotal  exchange  we  have: 


MINT(IPRIj  =  EPCI 
MEXT(EPCI)  =IPRI 
MINT(IPCI)  =  EPRI 
MEXT(EPRi)  =IPCI. 


We  cal!  this  U  pe  of  index  exchange  a  ‘'primary  exchange”.  For  abstrac¬ 
tion  purposes  we  define  a  function  Primary  Jndex_Exchange  (IPRI.IPCI)  which, 
given  the  current  tableau  representation  and  the  intrinsic  pivot  row  index  IPRI  and 
the  intrinsic  pivot  column  index  IPCl.  performs  the  appropriate  update  of  MINT() 
and  MEXT().  Ever\'  pivot  requires  a  primary  exchange. 

In  addition  to  the  partitioning  of  constraints  into  B  and  D.  we  are  re¬ 
quired  to  maintain  the  factorizations  ((?),(?j).  etc.)  within  each.  Having  performed 
the  primary  exchange,  it  is  possible  that  additional  exchanges  within  B  and/or  D 
are  required  to  restore  the  proper  factorization.  For  example,  if  the  primary  ex¬ 
change  remoN'e.-  an  /■’-tyf>e  constraint  from  D  (region  (t?))  and  places  it  in  region  (j) 
of  B.  an  additional  exrlianee  is  necessary  to  move  the  row  from  region  (j)  to  region 
(jj  I.  We  rei'er  to  smii  acticns  as  '■secondary  exchanges”.  Secondari'  exchanges  in  D 
corresjiond  to  row  exchanges  in  id.4)  and  secondary  exchanges  iti  B  to  column  ex- 
rhai.c  -  in  {■>.]  c  ']  hu-.  we  deline  the  function-^  Row  Jndex -Exchange!  IRll  .IIU'2  i 
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and  Column  Jndex_Exchange(ICIl,ICI2)  which  exchange  intrinsic  rows  IRIl  and 

IRI2  or  intrinsic  columns  ICIl  and  ICI2,  respectively.  Since  the  indexing  of  the 

RIN-O  array  corresponds  to  the  tableau  indexing,  we  assume  that  any  secondary 

exchanges  performed  on  the  tableau  are  also  performed  on  RIM(). 

Finally,  our  factorization  requires  that  the  factored  kernel,  Fu,  remain 

nonsingular.  It  is  possible  that  the  primary  and  secondary  exchanges  will  destroy 

the  nonsingularily  of  Fu-  When  this  occurs,  additional  row  exchanges  between 

regions  (i)  and  (Hi)  of  D  will  be  necessary  to  restore  the  required  nonsingularity  of 

Fu-  We  call  such  exchanges  “tertiary  exchanges”.  Again,  we  assume  that  tertiary 

exchanges  performed  on  the  tableau  are  also  performed  on  RIM(). 

The  factorizations  of  B  and  D  are  dynamic  in  nature,  meaning  that 

a  particular  pivot  may  cause  the  dimension  of  one  or  more  regions  of  B  and/or 

D  to  increase  or  decrease  by  one  row.  Such  dimension  changes  are  handled  by 

simply  adjusting  the  values  of  the  factorization  pointers.  We  define  the  functions 

Increment(XXX )  and  Decrement(XXX).  where  “XXX  is  MKC,  MFR,  MXR. 

NXR  or  NFR.  to  increment  or  decrement,  respectively,  the  appropriate  pointer 

value.  ,  ,  ,  , 

In  summary,  the  tableau  data  structures  are: 

MINT(] 

MEXT() 

R1M() 

ROWCOI-O 

.MKC.  MFR.  MXR. NFR.  NXR. 
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and  the  necessary  update  operations  are: 


Generate  JR,ow(7P/?/) 
Generate_Column(7PC/) 

Primary  Jndex_Exchange(7P7i7, 1  PCI) 
Row_Index_Exchange(77?71 , 77272) 

Column  Jndex_Exchange(7C7l,  7C72) 

Increment  ( A' XJV),  and 
Decrement(A' AA'). 

2.  The  Explicit  Transformation  Kernel 

As  can  be  seen  from  (5.4).  the  rows  of  correspond  to  region  (iit)  of 
the  tableau  (nonbasic  nonkey  variables)  and  the  columns  to  region  (J)  (basic  explicit 
constraints).  The  dimension  of  may  vary  as  the  algorithm  progresses,  and  the 
mechanism  for  these  changes  is  the  addition  and  deletion  of  rows  and  columns. 
W'c  thus  require  data  structures  that  permit  the  eflRcient  insertion,  deletion  and 
exchange  of  rows  and  columns  of  AJi  .  Further,  the  number  of  nonzero  elements 
in  varies  from  pivot  to  pivot,  independently  of  dimension  changes,  \^’e  store 
the.'C  nonzero  elements  in  a  stack,  referring  to  them  via  pointers  which  are  stored 
as  an  orthogonally  linked  list,  doubly  linked  by  row  and  doubly  linked  by  column. 
This  allows  the  efficient  update  of  and  accommodates  the  dynamics,  at  some 
expense  in  storage.  Contrary  to  the  tableau  representation,  we  maintain  the  explicit 
transformation  kernel  ref)resentat ion  indexed  extrinsically  rather  than  intrinsically. 
Thi''  allows  [iermiilat  ions  within  region  [j  ]  and/or  within  region  (??t)  of  the  tableau 


without  requiring  corresponding  permutations  in  the  columns  and/or  rows  of  the 
representation. 

We  again  define  a  suite  of  functions  which  operate  on  our  representation  of 
and  describe  its  update.  Add_Row(IRE)  and  Add_Column(IJE)  append  ex¬ 
trinsic  row  IRE  and  extrinsic  column  IJE,  respectively,  to  our  representation  of  Af/. 
We  will  ensure  that  the  current  representation  of  each  is  in  work  array  ROWCOL() 
prior  to  executing  these  functions.  Columns  will  always  be  appended  to  the  right- 
hand  side  of  Aj”/  and  rows  will  be  appended  to  the  bottom.  Delete  JRow(IRE)  and 
Delete_Column(IJE)  delete  extrinsic  row  IRE  and  IJE,  respectively,  from  the  rep- 
resentatio.i  of  Aj“/.  Replace_Row(IREl,IRE2)  causes  the  overlay  of  extrinsic  row 
IREl  by  extrinsic  row  IRE2.  The  representation  of  IRF  .  will  have  been  previously 
placed  in  ROWCOL().  Finallj’.  Upd?.te_Explic!t  TVansformation_Kernel  per¬ 
forms  a  pivot  update  of  the  representation  of  A^/,  using  the  current  representation 
of  the  pivot  row  and  pivot  column  in  ROWCOL(). 

3.  The  Factored  Kernel 

The  data  structures  and  update  actions  for  the  factored  kernel  vary  ac¬ 
cording  to  the  factorization  and  their  discussion  will  be  deferred.  Here  we  define  the 
necessary  abstract  update  functions. 

Is_Factored_KerneLSingular  tests  the  current  representation  of  Fn 
for  singularity  and  returns  the  appropriate  Boolean  value.  It  will  be  used  in  certain 
pivot  cases  to  determine  which  of  several  courses  of  action  should  be  taken. 

All  tertiary  exchanges  take  the  form  of  row  exchanges  between  regions 
(()  and  {lii}  of  (o.J).  Note  that  region  (?)  consists  of  the  nonnegativity  constraints 
of  key  variable-.  By  our  alternate  interpretation  of  the  nonnegativity  constraint 
indices  we  may  consider  these  tcj  be  key  variable  indices.  Since  the  columns  of  Fn 
are  ke\'  \'ariatile-.  wr 


interpret  region  (?)  as  containing  the  indices  of  the  columns 


of  /’ll.  An  exchange  of  rows  between  (i)  and  {Hi)  of  (5.4)  is  then  interpreted  as  an 
exchange  of  columns  of  Fn. 

Before  such  a  column  exchange  for  Fn  (and  thus  a  row  exchange  in  (5.4)) 
can  be  made,  we  must  frequently  identify  the  index  of  the  column  to  be  removed  from 
Fn  (which  is  an  index  in  region  (i)  of  (5.4))  and  the  index  of  the  column  which  will 
replace  it  (which  is  an  index  in  region  {Hi)  of  (5.4)).  We  thus  define  Find-Column 
_to_Remove(IOUT),  which  selects  from  among  the  indices  of  region  (i)  of  (5.4) 
the  intrinsic  index  lOUT  of  the  column  to  be  removed  from  Fn-  Similarly,  Find 
_Column_to_Add(IIN)  selects  from  among  the  indices  of  region  {Hi)  of  (5.4)  the 
intrinsic  index  IIN  of  the  column  to  be  be  added  to  F]].  Finally.  Update_Factored 
-Kernel  updates  the  data  structures  used  to  represent  Fji. 

E.  COMPLETE  ALGORITHM  DESCRIPTION 

We  are  now  in  a  position  to  fully  describe  our  implementation  of  the 
algorithm.  We  do  so  by  expanding  the  discussion  of  each  step  given  previously  in 
section  C. 

1.  Initialise  the  algorithm  using  the  origin  as  the  initial  basic  solution.  Then 

If 
If 

‘j.  Stop  if  tlu'  current  si)hition  is  optimal.  The  current  solution  is  optimal  if 

(RIMilHi  >  0  fur  \  <  IH  <  M  and  F{IM[JC)  <  0  for  .\/  +  l  <  JC  <  A/  +  .V. 

3,  (a)  Select  a  c  iulated  primal  oi  dual  constraint  as  the  target  row  or  column. 

respc’Ct  i\'el\ .  lor  purp()''C-s  of  exposition,  assume'  the  current  sedution  is 

''  i 


=  (-/) 


primal  feasible  and  we  are  executing  the  primal  algorithm.  Then  the 
target  row  is  the  bottom  row,  and  we  select  the  intrinsic  index  IPCl  of  a 
column  which  will  allow  us  to  make  a  gain  in  the  value  of  the  extremal 
function  (i.e.,  RIM{IPCI)  >  0). 

(b)  Generate_Column(IPCI). 

(c)  Calculate  ratios  by  computing  RIM{I R)l ROWCOL{I R)  for  those  1  < 
!R  <  M  with  ROWCOL{IR)  >  0.  M  R0WC0L{1  R)  <  0  for  all  7i?,  1  < 
1 R  <  A/  stop,  since  the  problem  is  primal  unbounded.  Otherwise,  select 
IPRI  to  be  the  intrinsic  index  yielding  the  smallest  such  ratio.  Break 
ties  in  accordance  with  the  following  hierarchy;  region  {Hi)  first,  then 
regie."  (7?),  then  (?)  and  finally  region  (iv).  Within  a  region,  break  ties 
by  selecting  IPRI  to  be  the  first  such  index  encountered. 

(d)  Generate_Row(IPRI ). 

4.  If.  contrary  to  our  assumpt  ion  in  2. a.  the  current  solution  is  not  primal  feeisible. 
the  target  row  would  be  some  row  IPRI  of  (5.4).  1  <  IPRI  <  M .  rather  than 
the  bottom  row.  W'e  would  proceed  by  next  executing  Generate  Row(lPRl). 
If  RO\\'COL{.JC)  >  0  for  all  A/  +  1  <  JC  <  M  +  A’,  we  would  conclude  that 
the  problenj  is  primal  infeasible. 

5.  (a)  Primary _Index_Exchange(IPRI.IPCi). 

fb)  Pivot  update  RIM()  and  OPT  using  ROWCOL(). 

(c)  Perform  the  necessary  secondary  and  tertiary  exchanges,  as  shown  in  Ta¬ 
ble  2.  N'r  .<■  that  some  tertiary  update  actions  are  conditional,  depending 
upon  the  singularity  of  Fj).  We  use  a  notation  similar  to  the  FORTR.-\.\ 
BLOCK  II  statement  to  indicate  the  conditional  actions. 


(d)  Update_Factored_Kernel. 

(e)  Update JExplicit -Transformation -Kernel. 

(f)  Go  to  Step  2. 

We  now  give  a  more  detailed  explanation  of  the  secondary  and  tertiary 
exchanges  listed  in  Table  5.2,  discussing  each  pivot  coordinate  individually. 

1.  Pivotal  coordinate  ((i),  (j)}.  EPCl,  the  extrinsic  index  of  an  E-type  constraint 
located  in  position  IPCI  of  £?,  is  exchanged  with  EPRI,  the  extrinsic  index 
of  a  nonnegativity  constraint  located  in  position  IPRI  of  D.  The  initial  pivot 
exchange  places  EPCI  in  position  IPRI  of  (2)  and  EPRI  in  position  IPCI  of 
(j),  and  thus  secondary  row  and  column  exchanges  are  necessary  to  restore  the 
structure  of  the  tableau.  Beginning  with  the  column  structure,  we  exchange 
column  EPRI  in  position  IPCI  of  region  [j]  with  the  extrinsic  index  located  in 
position  NXR  of  region  {jjj)-  and  then  decrement  the  value  of  NXR.  This  has 
the  effect  of  shifting  the  starting  column  index  of  region  (;;)  one  position  to 
the  left.  Now  extrinsic  index  EPRI  i.s  in  position  (NXR-!-])  of  the  tableau  and 
thus  is  still  misi)laccd.  Wo  then  exchange  EPRI  in  position  (NXR-i-1)  with  the 
extrinsic  index  located  in  position  NFR.  1  his  places  EPRI  contiguously  with 
region  {jjj)-  so  _  decrementing  NFR  we  restore  the  column  factorization. 
The  effect  has  bt'cn  to  decrease  the  dimension  of  (; )  by  one.  increase  the 
dimension  of  (jjj)  by  one  and  to  shift  the  entire  {jj)  region  one  position 
to  tlie  left.  In  the  row  structure.  EPCI  in  position  IPRI  of  region  (f)  must 
eventually  l.o  mo\ed  10  regioii  (rr).  However,  prior  to  this  pivot  EPRI  was  a 
key  variable.  IcKated  in  rc'gion  (?).  Its  removal  has  disturbed  the  structure  and 
t  hiis  t  he  rionsingulariiy  o  w  e  must  tlierefore  identih  a  column  currently 

in  reciof)  (///i  \^hieh  can  be  used  to  restore  the  nonsincularity  of  /  ]].  We 


thus  invoke  Find_Column_to_Add(IIN),  which  identifies  the  intrinsic  index 
IIN  of  such  a  column.  Suppose  the  extrinsic  index  of  the  column  in  position 
IIN  is  KX.  We  exchange  EPCI  in  position  IPRI  of  region  (i)  with  KX  in 
position  IIN  of  region  (ni).  We  now  have  extrinsic  index  EPCI  in  position 
IIN  of  {Hi)  and  we  must  move  it  to  region  (iv).  We  thus  exchange  EPCI 
in  position  IIN  in  region  {Hi)  with  the  extrinsic  index  in  position  NXR  of 
{Hi),  and  then  decrement  NXR.  This  completes  the  restructuring  ae  row- 
factorization  with  the  effect  of  decreasing  the  dimension  of  region  {Hi)  while 
increasing  that  of  region  {iv).  The  effect  of  these  exchanges  has  been  to  reduce 
the  dimension  of  through  the  deletion  of  column  EPCI  and  row  KX.  To 
maintain  the  proper  structure  of  v.e  therefore  Delete_Row(KX)  and 
Delete_Column(EPCl). 

2.  Pivotal  coordinate  ((;)•  {jj))-  EPCI,  the  extrinsic  index  of  a  basic  F-type  con¬ 
straint  located  in  position  IPCI  of  B.  is  exchanged  with  EPRI,  the  extrinsic  in¬ 
dex  of  a  nonbasic  nonnegativity  constraint  located  in  position  IPRI  of  F.  Since 
both  EPCI  and  EPRI  are  misplaced  after  the  primary  exchange,  secondary  row 
and  column  exchanges  are  needed.  In  the  columns,  we  exchange  EPRI  in  po¬ 
sition  IPCI  of  region  (j  ?)  with  ihe  extrinsic  index  located  in  position  NFR  of 
region  {jjl  Decrementing  NFR  then  completes  the  column  exchang-  s,  leaving 
the  dimen>ion  of  region  {j j  )  reduced  by  one  and  that  of  {jjj)  increased  by 
one.  In  the  rows,  we  exchange  EPCI  in  position  IPRI  of  region  (?)  with  the  ex¬ 
trinsic  indi’x  in  position  MKC  of  region  (?).  Decrementing  MKC  then  restores 
the  row  farton  d  .>-trurinie  of  the  tableau.  R  remains  to  be  seen,  however, 
whether  wha:  remains  of  / after  the  removal  of  row  EPCI  and  column  EPRI 
is  still  noii'-insjuhi;.  We  test  by  invoking  Is_Factored_KerneLSing-ular. 
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If  the  response  is  FALSE,  then  no  tertiary  exchanges  are  required.  If,  how¬ 
ever,  the  response  is  TRUE,  we  must  identify  a  column  exchange  that  we  will 
restore  the  nonsingularity  of  F^.  To  identify  the  two  columns  which  will  be  in¬ 
volved  in  the  exchange,  we  invoke  Find_Column_to_Reniove(IOUT)  which 
returns  the  intrinsic  index  lOUT  (with  associated  extrinsic  index  NX,  say)  of 
a  column  to  be  removed  from  Fn.  We  then  call  Find_Column_to_Add(IIN) 
which  identifies  the  intrinsic  index  IIN  (with  associated  extrinsic  index  KX, 
say)  of  a  column  in  region  (Hi)  which  can  be  used  to  replace  lOUT.  We  now  in¬ 
voke  Row  Jndex_Exchange(IIN,IOUT),  which  completes  the  update  of  the 
row  factorization.  This  row  exchange  results  in  the  replacement  of  row  KX  of 
by  row  NX  of  the  tableau.  We  currently  have  no  representation  of  row  KX  of 
the  tableau,  so  we  compute  it  by  invoking  Generate JRow(IOUT).  We  may 
then  restore  the  structure  of  AJ")'  by  invoking  Replace_Row(KX.NX). 

3.  Pivotal  coordinate  ((?)-(ijj))-  EPCl.  the  extrinsic  index  of  a  basic  nonneg¬ 
ativity  constraint  located  in  position  IPCl  of  B,  is  exchanged  with  extrinsic 
index  EPRI  located  in  position  IPRl  of  D.  The  primary  exchange  properly 
places  EPHl  in  region  {jjj).  and  thus  no  secondary  column  exchanges  are 
needed.  Likewise.  EPCl  is  properly  [)laced  in  position  IPRl  of  region  (?),  so 
no  secondary  row  exchanges  are  required.  Prior  to  the  pivot.  EPRI  was  desig¬ 
nated  a  key  column  and  its  removal  disturbed  the  structure  of  Fn.  Since  EPCl 
is  itself  a  column,  it  is  possible  that  it  may  be  a  replacement  for  EPRI.  To 
determine  if  thi^  is  thn  case,  we  invoke  Is_Factored_Kernel_Singular.  If  the 
respon^f  i^  E.M.SE.  tin  direct  exchange  of  column  EPCl  for  column  EPRI  in 
i*-  i>i'-tifi(  fl  a;.i!  t  lie  tatileau  update  is  complete.  If  the  response  is  TRUE,  how- 
e\-er.  a  lertiarx  lew  exchange  is  required  EPCl.  currently  located  in  pc)sition 


IPRI  of  region  (i),  will  be  the  column  we  remove  from  Fn,  and  its  replacement 
is  found  by  invoking  Find_Column_to_Add(IIN)  which  returns  the  intrinsic 
index  IIN  (associated  with  extrinsic  index  KX,  say)  located  in  region  (iii).  We 
perform  a  tertiary  exchange  by  invoking  Row J[ndexJExchange(IIN, IPRI). 
This  completes  the  row  factorization  update.  This  tertiary  exchange  results 
in  the  replacement  of  row  KX  of  by  row  EPCI  of  the  tableau.  We  must 
therefore  perform  the  corresponding  row  exchange  in  also.  To  perform 
such  an  exchange,  we  require  the  current  representation  of  row  IPRI  of  the 
tableau.  This  is  precisely  the  pivot  row  and  is  already  available  in  the  work 
array  ROWCOL().  Thus,  we  may  invoke  Replace_Row(KX,EPCl)  directly. 

Pivotal  coordinate  ((ii),(j))-  EPCi,  the  extrinsic  index  o^  a  basic  F-typc 
constraint  located  in  position  IPCI  of  B,  is  exchanged  with  EPRl,  the  ex¬ 
trinsic  index  of  a  nonbasic  F-type  constraint  located  in  position  IPRI  of  D. 
After  the  primary  exchange,  both  EPCI  (in  (i?))  and  EPRI  (in  (j))  are  mis¬ 
placed,  and  thus  both  secondary  row  and  column  exchanges  are  necessary. 
We  exchange  EPRI  in  position  IPCI  of  region  {j)  with  the  extrinsic  index 
in  position  NXR  of  region  (j)  and  then  decrement  NXR.  This  has  the  ef¬ 
fect  of  reducing  the  dimension  of  region  (j)  while  increasing  that  of  region 
(JJ).  Notice  that  we  have  added  row  EPRI  to  the  structure  of  Fn,  and 
thus  an  additional  column  must  at  some  point  be  identified.  Several  row 
exchanges  are  needed.  First,  we  exchange  row  EPCI  in  position  IPRI  of 
region  (n)  witli  the  extrinsic  index  located  in  position  (MKC+1)  of  region 
(2?)  and  then  decrement  MKC.  This  restores  the  structure  of  region  (n)  by 
reducing  iis  dimension  by  one  row.  We  now  must  execute  a  series  of  ex¬ 
changes  that  ultimatelv  places  EPCI  in  region  (u')  and  also  adds  a  column 


to  region  (t)  to  restore  the  nonsingularity  of  Fu-  To  do  this,  we  first  iden¬ 
tify  the  column  to  be  added  to  Fu  by  invoking  Find_Column_to_Add(IIN), 
which  identifies  the  intrinsic  index  IIN  in  region  (Hi)  (associated  with  an 
extrinsic  index  KX,  say)  which  may  be  added  to  Fu-  Notice  that  the  dimen¬ 
sion  of  Fii  is  increased  by  one  through  the  addition  of  row  EPRI  and  col¬ 
umn  KX.  We  then  perform  Row  JndexJExchange(IIN,MKC),  which  places 
KX  in  position  MKC  of  region  (i)  (correctly)  and  places  EPCI  in  position 
IIN  of  region  (in)  (incorrectly).  To  complete  the  tertiary  row  exchanges,  we 
invoke  Row Jndex_Exchange(IIN,MXR),  followed  by  Decrement(MXR). 
This  properly  places  EPCI  in  region  (in)  by  increasing  the  dimension  of  re¬ 
gion  (iv)  by  one  row  while  decreasing  that  of  region  (Hi).  Finally,  notice  that 
the  dimension  of  has  been  reduced  by  one  through  the  removal  of  column 
EPCI  and  row  KX.  To  update  the  representation  of  accordingly,  we  invoke 
Delete_Row’(KX)  followed  by  Delete.Column(EPCI). 

Pivotal  coordinate  ((H),  (jj))-  EPCI,  the  extrinsic  index  of  a  basic  /’-type 
constraint  located  in  position  IPCI  of  B,  is  exchanged  with  EPRI,  the  ex¬ 
trinsic  index  of  a  nonbasic  T-type  constraint  located  in  position  IPRI  of  D. 
The  primary  exchange  places  both  EPRI  and  EPCI  in  their  propei  positions, 
so  no  secondary  row  or  column  exchanges  are  needed.  However,  the  struc¬ 
ture  of  ha.s  been  altered  through  the  addition  of  row  EPRI  and  the  deletion 
of  row  EPCI.  To  determine  if  the  structure  and  nonsingularity  of  Fu  has  re¬ 
mained  intact,  we  invoke  Is_Factored_Kernel_Singular.  If  the  response  is 
F.'XLSL.  the  tableau  update  i.s  complete.  If  the  response  is  TRUE,  we  must 
identify  a  column  exchange  that  will  restore  nonsingularit\.  Thus,  we  invoke 
Find_CoIunin_to_Remove(10l  T)  which  returns  the  intrinsic  index  lOUT 


(associated  with  the  extrinsic  index  NX,  say)  of  a  column  located  in  region 
(t)  which  may  be  removed  from  Fn.  We  invoke  Find_Column_to_A.dd(IIN) 
which  returns  the  intrinsic  index  (associated  with  an  extrinsic  index  KX,  say) 
of  the  column  IIN  located  in  region  (tti)  which  may  be  added  to  Fu.  We  then 
perform  a  tertiary  exchange  by  invoking  Row  Jndex_Exchange(IIN,IOUT). 
This  last  action  results  in  the  replacement  of  row  IIN  of  by  row  lOUT  of 
the  tableau.  We  have  no  current  representation  of  row  lOUT,  and  to  compute 
one  we  invoke  Generate_Row(IOUT).  We  then  update  the  representation  of 
by  invoking  Replace -Row(KX, NX). 

6.  Pivotal  coordinate  ((n),  {jjj))-  EPCI,  the  extrinsic  index  of  a  basic  nonnega¬ 
tivity  constraint  located  in  position  IPCI  of  B,  is  exchanged  with  EPRI,  the 
extrinsic  index  of  a  nouud.sir  F-type  constraint  located  in  position  IPRI  of  D. 
The  primary  exchange  misplaces  both  EPRI  and  EPCI  and  thus  both  sec¬ 
ondary  row  and  column  exchanges  are  necessary.  We  exchange  EPRI.  located 
in  position  IPCI  of  region  (jjj)  with  the  extrinsic  index  located  in  position 
NFR.  By  then  incrementing  KFR,  we  restore  the  column  factorization  struc¬ 
ture  by  increasing  the  dimension  of  region  (jj)  while  reducing  that  of  region 
(jjj)-  Notice  that  the  addition  of  row  EPRI  to  the  structure  of  Fu  indicates 
that  an  a  corresponding  column  must  also  be  located.  To  restore  the  structure 
of  region  (ii).  we  exchange  EPCI  in  position  IPRI  of  region  (n)  with  the  extrin¬ 
sic  index  in  position  MXC  of  region  (f).  Incrementing  MKC  restores  the  struc¬ 
ture  of  region  (??)  by  reducing  its  dimension  by  one  row.  Since  EPCI  is  a  col¬ 
umn  which  we  have  now  placed  in  position  MKC,  it  is  possible  that  Fu  is  now 
nonsingular.  To  determine  this,  we  invoke  Is_Factored JKerneLSingular. 
If  the  response  is  I  .ALSL.  the  tableau  is  eomplete.  If.  on  the  other  hand. 


the  response  is  TRUE,  we  must  perform  a  tertiary  column  exchange  to  re¬ 
store  the  nonsingularity  of  Fn.  We  will  be  removing  EPCl,  currently  located 
in  position  MKC,  from  Fn  and  replacing  it  with  the  column  identified  by 
Find_Column_to_Add(IIN),  which  is  extrinsic  index  KX  located  in  posi¬ 
tion  IIN  of  region  (ni).  We  then  perform  Row  Jndex_Exchange(IIN,MKC), 
which  completes  the  tertiary  exchange.  This  exchange  implies  a  row  exchange 
in  also.  We  are  required  to  replace  row  KX  of  with  the  current 
representation  of  row  EPCI  of  the  tableau.  Since  IPRl  is  the  pivot  row,  the 
current  representation  of  EPCI  is  readily  available  in  ROWCOL().  We  may 
then  complete  the  row  exchange  by  invoking  Replace_Row(KX,EPCI). 

7.  Pivotal  coordinate  EPCI,  the  extrinsic  index  of  a  basic  E^type 

constraint  located  in  position  IPCl  of  B,  is  exchanged  with  EPRI,  the  extrin¬ 
sic  index  of  a  nonbasic  nonnegativity  constraint  located  in  position  IPRI  of 
D.  The  primary  exchange  misplaces  both  EPRI  and  EPCI,  and  thus  both 
secondary  row  and  column  exchanges  are  necessary.  We  exchange  EPRI.  in 
position  IPCI  of  region  (j).  with  the  extrinsic  index  in  position  NXR  of  re¬ 
gion  (j).  We  then  exchange  EPRI.  now  in  position  NXR  of  region  (;).  with 
the  extrinsic  index  located  in  position  NFR  of  region  We  complete  the 

column  update  by  decrementing  both  NXR  and  NFR,  which  has  the  effect  of 
reducing  the  dimension  of  region  {j)  by  one  column,  shifting  region  {jj)  one 
column  to  the  left  in  the  tableau  and  increasing  the  dimension  of  region  {jjj) 
by  one  column.  For  the  row  update,  we  first  exchange  EPCI,  in  p  on  IPRI 
of  region  (?;;).  with  the  extrinsic  index  located  in  position  MXR  of  region 
(ui).  By  tlien  decrementing  MXR  we  restore  the  row  factorization  by  increas¬ 
ing  the  dimension  of  region  (n-)  by  one  row  while  decreasing  that  of  region 


{Hi).  The  effect  on  has  been  to  reduce  its  dimension  by  one  through  the 
deletion  of  column  EPCI  and  row  EPRI.  The  corresponding  update  to  the 
representation  of  .4^/  requires  that  we  invoke  Deiete_Row(EPRI)  followed 
by  Delete_Column(EPCI). 

8.  Pivotal  coordinate  {(iii),(jj)).  EPCI,  the  extrinsic  index  of  a  basic  F-type 
constraint  located  in  position  IPCI  of  B,  is  exchanged  with  EPRI,  the  extrin¬ 
sic  index  of  a  nonbasic  nonnegativity  constraint  located  in  position  IPRI  of 
D.  The  primary  exchange  misplaces  both  EPRI  and  EPCI,  and  thus  both 
secondary  row  and  column  exchanges  are  necessary.  We  exchange  EPRI,  cur¬ 
rently  located  in  position  IPCI  of  region  (jj),  with  the  extrinsic  index  currently 
located  in  position  NFR  of  region  (jj).  By  then  decrementing  NFR,  w'e  re¬ 
store  the  column  factorization  by  decreasing  the  dimension  of  region  {jj)  while 
increasing  that  of  region  {jjj).  Note  that  the  removal  of  EPCI  from  the  struc¬ 
ture  of  Fii  implies  that  the  dimension  of  Fu  decreases  by  one,  and  thus  a 
corresponding  column  must  be  removed  from  Fu.  To  locate  this  column,  we 
invoke  Find_Column_to_Remove(IOUT),  which  returns  the  intrinsic  index 
lOUT  (associated  with  the  extrinsic  index  NX,  say)  of  a  column  located  in 
region  (?)  which  !iia\'  be  removed  from  Fn  while  allowing  the  remaining  struc¬ 
ture  to  be  noiisingular.  We  exchange  EPCI.  currently  located  in  position  IPRI 
of  region  (???).  with  NX.  currently  located  in  position  lOUT  of  region  (?).  This 
restores  the  row  structure  of  region  {Hi),  but  EPCI  is  misplaced  in  region  (?). 
Therefore,  we  exchange  EPCI  in  position  lOUT  of  region  (?)  with  the  extrin¬ 
sic  index  located  in  position  MKC  of  region  (?).  By  then  decrementing  MKC, 
we  restore  tin'  st'urture  of  region  (?)  by  decreasing  its  dimension  and  that 
of  region  (?-  !  by  increasing  its  dimension.  The  exchange  of  NX  in  IlN  with 


EPCI  in  IPRI  implies  a  row  exchange  in  .  To  perform  this  exchange  in 
the  representation  of  we  must  replace  row  EPRI  of  with  the  current 
representation  of  row  NX  of  the  tableau.  Since  this  current  representation  is 
not  available,  we  invoke  Generate JFlow(IIN)  to  compute  it.  We  then  update 
our  representation  of  by  invoking  Replace_Row(EPRI,NX). 

9.  Pivotal  coordinate  EPCI,  the  extrinsic  index  of  a  basit,  nonneg¬ 

ativity  constraint  currently  located  in  position  IPCI  of  B,  is  exchanged  with 
EPRI,  the  extrinsic  index  of  a  nonbasic  nonnegativity  constraint  currently  lo¬ 
cated  in  position  IPRI  of  D.  The  primary  exchange  properly  places  EPRI  and 
EPCI,  so  no  secondary  exchanges  are  necessary.  Since  is  unaffected,  no 
additional  action  is  required. 

10.  Pivotal  coordinate  EPCI.  the  extrinsic  index  of  a  basic  E-type 

constraint  currently  located  in  position  IPCI  of  B,  is  exchanged  with  EPRI,  the 
extrinsic  index  of  a  nonbasic  £'-type  constraint  currently  located  in  position 
IPRI  of  D.  The  primary  exchange  properly  places  EPRI  and  EPCI.  so  no 
secondary  exchanges  are  necessary.  Fu  is  unaffected,  so  no  additional  action 
is  required. 

11.  Pivotal  coordinate  {{n-],{jj)).  EPCI.  the  extrinsic  index  of  a  basic  F-type 
constraint  currently  located  in  position  IPCI  of  B,  is  exchanged  with  EPRI. 
the  extrinsic  index  of  a  nonbasic  E-type  constraint  currently  located  in  po¬ 
sition  IPRI  of  D.  The  primary  exchange  misplaces  both  EPRI  and  EPCI. 
so  both  secondary  row  and  column  exchanges  are  necessary.  EPRI.  cur¬ 
rently  in  position  IPCI  of  region  (jj).  is  exchanged  with  the  extrinsic  in¬ 
dex  located  at  [losition  (N.XR  +  ll  of  region  {jj).  The  column  structure  is 
restored  by  incrementing  NXR.  increasing  the  dimension  of  region  (j)  while 


decreasing  that  of  region  {jj).  Note  that  since  the  dimension  of  region  (j) 
increases,  the  dimension  of  will  increase  also.  Also,  the  removal  of  EPCI  from 
-^11  implies  that  the  dimension  of  Fn  will  decrease.  Thus,  a  column  of  Fji 
must  be  identified  for  removal  as  well.  To  find  such  a  column,  we  invoke 
Find_Column_to_Remove(IOUT),  which  returns  the  intrinsic  index  lOUT 
(associated  with  the  extrinsic  index  NX,  say)  of  a  column  whose  removal  from 
Fii  allows  the  remaining  structure  of  Fu  to  be  nonsingular.  We  perform  a  ter¬ 
tiary  exchange  by  exchanging  NX,  currently  located  in  position  lOUT  of  region 
(f),  with  the  extrinsic  index  currently  located  in  position  MKC  of  region  (i). 
The  structure  of  region  (i)  is  then  restored  by  decrementing  MKC,  reducing 
the  dimension  of  region  (?)  while  increasing  that  of  region  (n).  To  properly  po¬ 
sition  EPCI  and  restore  the  structure  of  region  (n),  we  exchange  NX,  currently 
located  in  position  (MKC+1)  of  region  (??),  with  EPCI,  currently  located  in 
position  IPRI  of  region  (iv).  This  leaves  only  NX  misplaced.  We  thus  exchange 
NX.  currently  located  in  position  IPRI  of  region  {iv),  with  the  extrinsic  index 
located  in  position  MXR  of  region  (??).  The  row  update  is  completed  by  incre¬ 
menting  MXR,  increasing  the  dimension  of  region  {Hi)  while  decreasing  that 
of  region  {iv).  Since  the  dimension  of  Aj'/  has  been  increased,  we  must  up¬ 
date  our  representation.  We  are  required  to  add  the  current  representation  of 
column  EPRl.  which  is  available  as  the  pivot  column  in  ROWCOL().  We  also 
require  the  current  representation  of  row  NX  of  the  tableau,  which  is  not  cur¬ 
rently  available.  We  thus  invoke  Generate_Row(IOUT)  W'ith  these  tw’o  rep¬ 
resentations,  we  may  then  invoke  Add-Row(NX)  and  Add_Column(EPCl) 
to  complete  the  update  of  the  repre.sentation  of  Aj"/- 


12.  Pivotal  coordinate  {{iv).,{jjj))-  EPCI,  the  extrinsic  index  of  a  basic  nonneg¬ 
ativity  constraint  currently  located  in  position  IPCI  of  B,  is  exchanged  with 
EPRI,  the  extrinsic  index  of  a  nonbasic  E-type  constraint  currently  located  in 
position  IPRI  of  D.  The  primary  exchajige  misplaces  both  EPRI  and  EPCI. 
To  restore  the  column  factorization,  we  first  exchange  EPRI,  located  in  posi¬ 
tion  IPCI  of  region  (j  jj),  with  the  extrinsic  index  located  in  position  (NFR-fl) 
of  region  [jjj)-  We  then  exchange  EPCI,  now  in  position  (NFR+l)  of  region 
ijjj)i  with  the  extrinsic  index  located  in  position  (NXR+l).  We  complete 
the  column  update  by  incrementing  both  NXR  and  NFR,  which  has  the  effect 
of  increasing  the  dimension  of  region  (j),  shifting  region  (jj)  one  position  to 
the  right  and  decreasing  the  dimension  of  region  (jjj).  Since  the  dimension  of 
region  (j)  has  been  increased,  the  dimension  of  will  increase  as  well.  To 
restore  the  row  factorization,  we  exchange  EPCI,  currently  in  position  IPRI 
of  region  (zc),  with  the  extrinsic  index  located  in  position  (MXR+l)  of  re¬ 
gion  (rc).  The  row  update  is  completed  by  incrementing  MXR,  which  has 
the  effect  of  increasing  the  dimension  of  region  (n’t)  while  decreasing  that  of 
region  (iv).  The  dimension  of  has  been  increased  through  the  addition  of 
column  EPCI  and  row  EPRI.  Since  IPCI  is  the  pivot  column  and  IPRI  is  the 
pivot  row.  eacli  is  already  available  in  ROWCOL().  Thus  we  simply  invoke 
Add_Row(EPRI)  followed  by  Add_Column(EPCI)  to  complete  our  update 
of  the  representation  of  .dj”/. 


VI.  FACTORIZATION  OF  GENERALIZED 
UPPER  BOUND  ROWS 

A.  INTRODUCTION 

We  now  develop  the  first  of  three  speciahzations  of  the  factorization  approach. 
We  are  still  interested  in  the  problem: 

(GUB)  min  :  wy 

s.t.  :  Fy  <  b  F  p  by  n 

Ey  <  r  E  m  by  n 
—  7j/  <  0  —  7  n  by  n 

where  F,  F,  —7.  ic,  b  and  r  are  as  before.  We  now  require  that  the  F-type  constraints 
are  generalized  upper  bound  (GUB)  constraints.  Define  S,.i  =  1 _ ,p  to  be  pair- 

P 

wise  disjoint  subsets  of  the  set  A'  =  {1 . n}  and  further  define  Sc,  =  A'  -  (j  5, 

«=i 

(5c  may  be  empty).  Then  5.05,-  =  0  for  0  <  i  <  p.  0  <  ?'  <  p.  i  ?'  and 

p  ^ 

(J  5,  =  A'.  Gl'B  constraints  are  of  the  form 

1=0 


Y,  ^  1 . P-  K  =  ±1  (6-1) 

The  sets  5,.  ?  =  0 . p  are  called  GUB  sets  and  ?  is  the  GUB  set  index.  Sc 

identifies  tliose  variables  that  have  no  coefficient  in  any  GUB  constraint.  The  F- 
type  constraints  are  of  arbitrary  form. 

A  specialization  of  the  Simplex  Method  to  handle  GUB  constraints  was  first 
developed  by  Dantzig  and  \'an  Slyke  [19f)7j.  They  introduce  a  specialization  of 
(GIB)  with  5,  =  ]  and  strict  equality  constraints.  Their  algorithm  requires  a 
working  basis  of  dimension  (n?  4-  1).  McBride  [1972]  develops  a  specialization  of 
the  mutual  primal-dua!  method  to  solve  a  second  variant  of  (GUB)  with  5;  =  il 
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and  strict  equality  constraints.  A  dynamic  working  baisis  whose  dimension  cannot 
exceed  m  is  required.  A  computational  analysis  by  McBride  [1972]  predicts  that  the 
performance  of  this  algorithm  should  exceed  that  of  the  Simplex  Method  when  the 
proportion  of  GUB  constraints  in  the  m.odel  exceeds  approximately  29%,  and  his 
empirical  results  tend  to  confirm  this  analysis. 

VVe  extend  this  work  by  developing  a  specialization  of  the  factorization  ap¬ 
proach  which  does  not  require  all  GUB  constraints  to  be  binding.  We  will  see  that 
there  are  computational  efficiencies  to  be  gained  by  this  approach.  We  first  present 
the  factored  tableau  for  (GUB).  We  then  discuss  the  important  computational  is¬ 
sues. 


B.  THE  FACTORED  TABLEAU 

Recall  from  Chapter  -1  the  primal  row  basis  at  an  arbitrary  point  in  the  solution 
process  is: 


k 

/ 

-- 

n  —  i—k 

\ 

(./ ' 

i  1 1 

F,: 

Eyj 

Fn 

^  13 

(.7.77) 

\  6 

0 

) 

}n  -  I  -  k- 


(6.2) 


We  saw  in  Chapter  1  that  it  is  always  possible  to  identify  among  the  columns  of 
B  a  set  of  k  column^  such  that  the  submatrix  Fn  is  nonsingular.  Since  the  row,'- 
of  i/"!!  B]2  ■f'l''  corre--p<)iic!  to  GUB  constraints,  their  form  is  similar  to  (6.1).  with 
column  permut  a:  i' >ii-  ac(  uuntmg  for  the  difference.  Each  column  of  [Fn  Fij  F12]  is 
cither  zero,  or  a  sima'h  urn'  vecou.  and  thus  by  row  permutations  of  these  F-type 
constraints,  /  n  rna;.  In-  plaenl  into  the  form  of  a  k  by  k  signed  identity  matrix  An. 


U(. 


d  he  resultinc  KU  B'  row  basi'-  i>; 


k 


n— /— it 


B  = 


O') 

Oj) 

00) 


/ 

Ei\  £12 
•^11  Bi2 

V  0  0 


\ 

£13 

/'l3 

-I ) 


)l 

}n  —  I  —  k 


(6.3) 


The  corresponding  nonbasic  constraint  rows  are  then: 


D  = 


(0 

/ 

-1 

0 

“  1 

)k 

(??) 

0 

F 22 

F 23 

]p-k 

(???) 

0 

-/ 

0 

(?r) 

F21 

F22 

.^23  j 

}  m  —  1 

(6.4) 


Recall  that  wc  wish  to  observe  a  one-to-one  association  between  the  variables 
whoso  nonnegativity  constraints  (which  appear  in  region  (?)  of  D)  are  nonbasic 
(nonbinding)  and  the  F-typc  constraints  in  region  [jj)  of  B.  This  association  is 
exactl}'  analogous  to  the  idea  of  ‘'key"  variables  proposed  by  Dantzig  and  Van  Slyke 
[1967].  Wo  call  tho  variables  whoso  nonnegativity  constraints  appear  in  region  (?) 
of  D  "b'y"  variables,  and  there  is  one  such  variable  for  each  currently  basic  F- 
type  constraint  appearing  in  region  (_;_;)  of  B.  The  variables  whose  nonnegativity 
constraints  appear  in  region  (???)  of  D  are  then  called  “non-key". 

Tsing  (6.3)  and  ((i.  l)  and  recalling  the  expression  for  the  explicit  transforma¬ 
tion  kernel; 


-111  ”  ~  .^11^11  ^12)  —  (/'i:  —  TiiA]]/ 


the  principal  part  of  the  factored  tableau  is: 


Fxi  “■ 

AJ'j't  £13  “  £n*^u£is) 
u  -  £n.iu£i3  “ 

+  £n‘i»i£’ii'^ri‘t£i3  -  £u^it£is) 


C.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 
EXPRESSIONS  IN  COLUMN  GENERATION 

'Jo  exanniic  the  actions  required  by  Generate-Column(LC).  suppose  we 
want  to  generate  roluinn  LC  of  (he  principal  part  of  the  tableau  (6.5)  and  place 
the  result.'^  in  the  \'e('tor  =  (-!•  -z-  -a-  which  is  partitioned  to  conform  to 

(6.5),  Itewriting  (6.5|  in  a  more  convenient  manner,  we  have; 


ij) 

(■) 

r 

—  An  F 11  Mu' ] 

-  AnTij  £^11*^111  + 

(h) 

-Fn  l-^u  j 

~Fj] 

(...) 

An 

(.0 

-En 

IaJ",'!  -  £31  [“ *^11  ^nMii  jj 

-Ejj  [-AfiTiiAiij 

"■Eol  ^ 


-AiiFii  |.4n  (£i3  -  Ex\^\\Fij)\  +  Ana'll 

—  Fi2  [A^i'lE^ts  —  ^iiAiiFij)]  +  Fji 
AI',’{£ij  —  EiiAiiru' 

—Bji  I'^r/t^is  "  £\x^\\F iji] 

—  |— Aufdi  A7i'(i^i3  “  Exx-^mF 13))  ■*■  ^11  ^"11]  ■*■  Fiz  } 


(6.6) 

To  highliglit  tlio.'^e  computations  involving  tcrtns  consisting  cntirf’ly  of  zeros  and  plus 
and  minus  ones,  we  introduce  the  following  notation:  a  general  matrix  product  is 
denoted  by  as  in  .•Ij','  ■  Tni-  A  ‘■sim{)le'’  matrix  product  (i.e.,  one  in  which  one 
of  (he  terms  consists  entirely  of  zero'-',  fdus  ones  and  minus  ones)  is  denoted  by  o,  as 
in  £11  o£i2.  1  he  calculations  may  then  proeex'd  as  follows  (note  that  is  a  unit 
vector  with  a  plus  one  in  |>usition  1,(.'): 


1 .  Calculate  region  ( lu  ): 

(a)  If  column  IJ'  i--  in  (/).  (hen: 


=  A,; 


o  c 


rr  ^  ^  A-i)rr 


Ml 


(b)  If  column  !,<  '  i-  in  (_/ ; ).  then; 


~ •  i 1 1'  ■  (  f-ii  '  (  A] 


<11 


(c)  If  column  LC  is  in  {jjj),  then: 


Z3  =  ■  [(£13)^^  -  E,,  o  A„  o  (Fj3)^^J 

Notice  that  is  either  null  or  it  is  a  signed  unit  vector,  say  in  row 

k,  in  which  case  En  o  An  o  =  ±(£n)^. 

2.  Calculate  region  (z): 

(a)  If  column  LC  is  in  (j),  then: 


Z\  —  —All  o  Fi2  o  Z3 

(b)  If  column  LC  is  in  (jJ).  then: 

Zj  =  —All  ^  ^12  *3  + 

(c)  If  column  LC  is  in  (jjj).  then: 


—  —All  o  Fi2  o  Z3  -i-  All  o  (Fis)^^. 


Notice  that  the  computation  —F12OZ3  involves  only  additions  and  sub¬ 
tractions.  and  that  each  case  in  region  (7)  differs  only  by  an  addition 
suffix. 

3.  Calculate  regicm  (;i  ): 


100 


(a)  If  column  LC  is  in  (j),  then: 


Zi  —  ~F22  ^  23 


(b)  If  column  LC  is  in  {jj),  then: 


Z2  =  -F22  O  23 


(c)  If  column  LC  is  in  {jjj),  then: 


22  —  —F22  o  23  +  (-^23)^^ 


4.  Calculate  region  (zr): 


(a)  If  column  LC  is  in  (j).  then: 


24  =  — 


21 


^22  •  ~3 


(b)  If  column  LC'  is  in  {jj).  then: 


24  —  —£^21  2\  —  E\ 


22  •  ‘3 


(c)  If  column  LC  is  in  (jjj).  then: 


'■)  —  “■£21  ■  2]  —  Lo;  •  23  -f  (£2-C 
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Then  column  LC  of  (6.6)  is  available  as  =  {zi,  22)  ^3,  Notice  that 

floating  point  multiplications  and/or  divisions  are  necessary  only  for  computation  of 
regions  (in)  and  (iu);  also  note  that  regions  (i),  (it)  and  (iv)  are  linear  combinations 
of  region  (iii)  and  one  another,  and  that  each  term  in  these  regions  differs  only  by 
an  additive  suffix. 

D.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 
EXPRESSIONS  IN  ROW  GENERATION 

Now  we  examine  the  actions  of  Generate_Row(LR).  Suppose  we  desire  to 
generate  row  LR  and  place  the  results  of  the  computation  in  the  vector  z  =  (is,  ie*  ^7) 
which  is  partitioned  to  conform  to  (6.5).  Rewriting  (6.5)  in  a  convenient  form: 


o  U'  (;;j} 


(•) 

j-Ajj£j).4J'j’l  £j3  i- 

|•-Al]£ £j|  /  /Alii  £ij 

(l.i 

"  {  ”  •  £u  ) 

1  £,3  *  |-  ■ 

[  £ij /  An  1  £13  *  £1; 

(ill) 

•ir; 

—  1 » A,,'  1  £j,  }  All 

j A,j’ '  £,5  *  i 

-(iAr/i£„]A„}£n 

(..I 

-  1 1  ~  £■::  £u  -*■  £j.  f  Ai 

1  |f£jiAn 

£lJ  ”  £33  *A,;'  £13 

(6.7) 


IIL- 


Then  we  may  proceed  as  follows  (note  that  em  is  a  unit  vector  with  a  plus  one  in 
position  LR): 

1.  Calculate  region  (j): 

(a)  If  row  LR  is  in  (t),  then: 

h  =  [—(^11  )lr  o  -fiz]  o 

(b)  If  row  LR  is  in  (it),  then: 


is  =  (  — ^22)Lft  O 

(c)  If  row  LR  is  in  (fn).  then: 

is  =  =  (^n 

(d)  If  row  LR  is  in  (u-),  then: 

■^5  =  [(■£'21  )lr  o  -^11  —  {£'22)lr]  •  •^ii' 

Notice  that  the  computation  of  the  term  (.£21  )l/?  ^  -^ii  ^  ~  (£22)1;? 

involves  only  additions  and  subtractions. 

2.  Calculate  region  [jj]: 

(a)  If  row  LR  is  in  (?).  then: 


•^1.  =  I  — •  £11  +  cl.v)  ^  -N)) 
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(b)  If  row  LR  is  in  (ii),  then: 


^6  —  (“25  •  ■£'11 )  O  ^11 


(c)  If  row  LR  is  in  (in),  then: 


26  —  (“-^S  ‘  ■£'ll)  O  ^11 


(d)  If  row  LR  IS  in  (in),  then: 


■^6  —  (“■^s  •  ■E'n  “  (•£'21  )l;?)  An 


3.  Calculate  region  Ujj)- 


(a)  If  row  LR  is  in  (i),  then: 


i?  =  25  •  E\3  zeO  Fi3 

(b)  If  row  LR  is  in  (n),  then: 


Z~  —  Z3  ■  E]3  -I-  Ze  O  E]3  -i-  {F33)ir 


(c)  If  row  LR  is  in  {Hi),  then: 


27  —  Z^  ■  E\3  Z(,  O  F\3 
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(d)  If  row  LR  is  in  (tv),  then: 


Zt  —  Zs  ■  E\3  +  26  O  Fi3  +  (-£^23)lR- 

The  row  LR  of  (6.7)  is  available  as  f  =  (25,  ie,  27).  Note  that  floating  point 
multiplications  and/or  divisions  are  avoided  for  terms  involving  Fu,  F22  and  F23, 
that  regions  (jj)  and  (jjj)  are  linear  combinations  of  region  (j),  and  that  many 
ca^es  differ  only  by  additive  infix  terms. 

This  computational  scheme  extends  the  approach  of  the  mutual  primal-dual 
method  by  introducing  additional  linear  dependencies  into  the  tableau  and  exploits 
those  dependencies  to  improve  efficiency.  Empirical  evidence  suggests  that  this 
approach  also  improves  numerical  stability  (McBride  [1989]). 

E.  DATA  STRUCTURES 

The  essential  information  contained  in  the  factored  kernel  when  the  factored 
rows  are  GUB  constraints  is  simply  the  unique  mapping  between  each  basic  factored 
constraint  and  its  corresponding  ‘‘key"  nonbasic  variable.  We  require  a  represen¬ 
tation  of  this  mapping  which  is  compact,  can  be  efficiently  accessed  and  is  easily 
updated. 

One  possibility  is  to  maintain  the  intrinsic  ordering  of  the  tableau  in  such  a  way 
that  the  factored  row/key  variable  relationship  can  be  derived  implicitly.  Region 
(jj)  of  the  tableau  corresponds  to  the  basic  factored  constraints  of  B.  Region  (?) 
of  the  tableau  correspond^  to  the  nonnegativity  constraints  associated  with  the  key 
variables.  If  we  maintain  the  ordering  of  [jj)  and  (?)  so  that  the  position  of  (?) 
contains  the  index  of  the  key  variable  associated  with  the  factored  constraint  whose 
index  is  in  the  pn-itron  of  (jj).  we  ha\-e  an  implicit  representation  of  Fij. 


Two  difficulties  arise  with  such  an  approach  which  lead  us  to  choose  an  alter¬ 
nate  representation.  First,  such  an  ordering  complicates  the  tableau  update  scheme 
presented  in  Chapter  5  and  potentially  increases  the  computational  expense  of  all 
tableau  updates  which  affect  Fu.  Second,  several  update  cases  require  the  iden¬ 
tification  of  GUB  set  membership  for  variables  which  are  not  currently  assigned 
to  region  (t)  and  thus  are  not  currently  part  of  the  Fu  structure.  This  may  be 
done  by  accessing  the  original  problem  data  by  column  and  scanning  that  column 
for  the  index  of  the  GUB  row,  if  any,  in  which  the  column  has  a  nonzero  element. 
Since  the  original  problem  data  is  stored  in  a  super-sparse  form,  this  scheme  re¬ 
quires  indirect  computer  memory  addressing.  We  prefer  a  method  which  requires 
less  computational  expense. 

We  thus  introduce  an  additional  data  structure  to  manage  the  structure  of 
Fii-  KEY(IJ)  is  an  array  of  type  INTEGER,  dimensioned  from  1  to  (m  -I-  n), 
which  for  each  basic  factored  constraint  records  the  extrinsic  index  of  the  associated 
key  variable,  and  for  each  nonbasic  key  variable  records  the  extrinsic  index  of  the 
corresponding  basic  factored  constraint.  Additionally,  for  each  variable  which  is  not 
key  (i.e.,  either  non-kpy  variables  or  btisic  variables),  KEY(IJ)  records  the  extrinsic 
index  of  the  factored  constraint  for  which  it  may  serve  as  a  key  variable.  If  we 
interpret  the  extrinsic  index  of  the  factored  constraints  as  GUB  set  indices,  then 
for  every  variable  KE^'(I>1)  contains  the  GUB  set  index.  Note  that  since  GUB 
set  membership  is  fixed  for  a  given  variable,  all  column  indices  KEY(IJ)  can  be 
initialized  once  and  need  never  be  updated.  Only  the  factored  constraint  indices 
require  updating. 

F.  FACTORED  KERNEL  UPDATE  ACTIONS 

Recall  from  Chapter  -5  the  details  of  the  factored  kernel  update  actions 
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Is_Factored_Kernel_Singular, 
Find_Column_To_Remove(7FA'C), 
Find_Column_To_Add(//’A'C),and 
Update-Factored -Kernel, 

are  factorization-specific.  We  now  discuss  these  details  for  the  GUB  factorization. 

The  row-permuted  diagonal  structure  of  Fn  (An)  imphes  that  the  question  of 
singularity  arising  from  rank-one  updates  can  be  resolved  strictly  on  the  basis  of  local 
information.  That  is,  the  exchange,  deletion  or  addition  of  a  row  leads  immediately 
to  the  identification  of  a  unique  GUB  set.  The  search  for  a  corresponding  column 
(variable)  to  complete  the  exchange,  deletion  or  addition  is  then  limited  to  those 
sharing  this  GUB  set  membership.  The  simplicity  of  this  structure  greatly  enhances 
the  efficiency  of  the  necessary  update  actions. 

For  example,  consider  pivot  update  case  ((n),(ii)),  in  which  basic  factored 
constraint  EPCI.  located  in  tableau  position  IPCI  is  removed  from  the  row  basis 
and  nonbasic  factored  constraint  EPRI,  located  in  tableau  position  IPRI  replaces  it. 
Because  EPCI  was  basic,  there  is  a  key  variable,  say  EPCIKV,  located  in  position 
IPCIK\'  of  region  (?)  of  the  tableau.  Because  GUB  sets  are,  by  definition,  pair¬ 
wise  disjoint,  EPCIK\'  may  not  serve  as  a  key  variable  for  the  new  basic  factored 
constraint  EPRI.  Thus,  a  replacement  variable  among  those  in  region  (Hi)  of  the 
tableau  must  be  found.  Such  a  replacement  must  exist,  for  otherwise  B  is  singular. 
To  locate  such  a  variable,  we  scan  the  variables  in  region  (n?),  searching  for  one 
which  is  a  member  of  GUB  set  EPRI.  The  test  is  simply  this:  for  each  variable  JC 
in  region  (???).  if 

KEY(JC)  =  EPRI. 


107 


then  accept  JC  as  the  key  variable  for  EPRI  and  proceed  with  the  secondary  tableau 
update  exchange.  Otherwise,  continue  the  scan. 

Thus,  the  action  Find_Column_To_Add(EKV)  scans  the  indices  of  region 
(in'),  performing  the  GUB  set  membership  test  on  eaeh  index.  The  first  such  index 
for  which  the  test  is  true  is  chosen  as  EKV. 

The  action  Find-Column _To_Remove(EK\^)  is  required  when  a  constraint 
ERI  is  removed  from  Fn-  EKV  is  then  the  key  variable  corresponding  to  ERI  and 
is  determined  by; 


EKV  =  KEY(ERI). 

Is_Factored_KerneLSingular  is  required  when  the  constraint  EPRI  is  to  be 
added  to  Fn  and  the  pivot  column  index  EPCI  corresponds  to  a  variable.  Thus, 
EPCl  is  a  possible  key  variable  for  EPRI  and  is  a  convenient  candidate.  It  may  be 
accepted  as  the  key  variable  for  EPRI  if  and  only  if  Fn  remains  nonsingular  after  the 
addition  cf  both  EPRI  and  EPCI,  and  this  is  the  case  if  EPCI  is  a  member  of  GUB 
set  EPRI.  Thus.  Is_Factored_KerneLSingular  simply  compares  KEY(EPCI)  and 
EPRI.  If  they  are  unequal,  Fu  is  singular  and  the  result  “true’'  is  returned.  Other¬ 
wise.  the  result  is  "false". 

Finally,  we  consider  Update_Factored_Kernel.  Since  the  GUB  set  mem¬ 
bership  of  a  variable  does  not  change,  only  constraint  index  updates  of  KEV(IR) 
are  required.  For  example,  if  the  old  key  variable  index  for  constraint  EPRI  is 
EPRIKVO  and  the  new  index  is  EPRIKV'N,  the  required  update  is: 

KEY(EPRI)  =  EPRIKVN. 
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VII.  FACTORIZATION  OF  PURE  NETWORK 

ROWS 


A.  INTRODUCTION 

We  remain  interested  in  the  problem: 

(PNSC)  min  :  wy 

s.t.  :  Fy  <  b  •  F  p  by  n 

Ey  <  r  ;  E  m  by  m 

—ly  <  0  ;  —I  n  by  m 

where  F,  E,  — /,  il\  b  and  r  are  as  before.  We  now  require  that  each  column 
of  F  consist  of  zero,  one  or  two  nonzero  elements  If  a  column  contains  a  single 
nonzero  element  within  the  rows  of  F,  it  may  be  either  a  plus  one  or  a  minus 
one.  If  a  column  contains  two  nonzero  elements  within  the  rows  of  F,  one  must 
be  a  plus  one  and  the  other  a  minus  one.  F-type  constraints  may  also  be  in¬ 
terpreted  as  defining  the  node-arc  incidence  matrix  of  a  directed  graph.  Define 

J  =  {0. . . . , p}  ,  =  {I . n } ,  and  A'  =  {n,.  ?  6  F  —  {0} }  to  be  a  set  of  nodes  and 

A  =  €  V.F  I  a*'  =  (tj,,  Uj).?  G  I.j  €  to  be  a  set  of  arcs  with  ordered  pairs 

of  nodes  (tail. head)  as  elements  indexed  by  k.  Note  that  we  interpret  node  0  as  a 

null  node,  and  thus  an  arc  incident  to  node  0  is  viewed  a-s  a  single-ended  arc.  Then 

a  graph  is  defined  a,*;  ^  =  {A  .  .4},  The  node-arc  incidence  matrix  of  is  a  matrix 
F  consisting  of  p  rows  (one  for  each  node  in  A  )  and  n  columns  (one  for  each  arc  in 
A)  with  elements: 


i-fl  if  (7^  =  (v,.rij)  for  sonic  cA 

—  1  if  o*  =  for  some  G  A' 

(I  otherwise 


Thus,  if  arc  a*  =  (n,,nj)  is  represented  by  column  k  of  F,  then 

0  \ 

■  row 

1 
0 

_  1  j  row 

0 

Oy 

and  each  column  of  F  is  either  all  zeros,  contains  exactly  one  nonzero  element  (which 
may  be  either  plus  one  or  minus  one)  or  contains  exactly  two  nonzero  elements  (a 
plus  one  and  a  minus  one). 

A  specialization  of  (PXSC)  that  has  been  widely  studied  arises  when  the  F~ 
type  constraints  are  taken  to  be  equalities  and  the  £-type  constraints  are  vacuous, 
resulting  in: 

(PNE)  min  :  U'lj 
y 

s.t .  :  Fy  =  h  F  p  by  ti 

£  0  —  I  n  by  n 

\ery  efficient  sj'ecializations  of  the  primal  Simplex  Method  have  been  developed 
and  implemented  to  solve  fPNE)  (e.g..  Srinivasan  and  Thompson  [197.3j.  Glovei.  et. 
al.  [(1974).  Bradley.  Brown  and  Graves  [1977]).  These  algorithms  exploit  some  well 
known  properties-  of  I  (we  assume  that  one  redundant  row  has  been  removed  from 

/•'I: 

1.  Every  primal  Simplex  Ijasis  B  of  F  (consisting  of  (p-  1)  linearly  independent 
column-  of  /■  i  can  f>e  triangulated  tw  row  and  column  permutation-. 


2.  Every  such  basis  B  is  itself  the  node-arc  incidence  matrix  of  a  subgraph  of  Q, 
and  this  subgraph  is  a  rooted  spanning  tree,  and 

3.  F  is  totally  unimodular,  implying  that  if  th^  and  65  are  integer  {p  —  l)-vectors 
and  Xi  and  X2  are  {p  —  1)- vectors  of  unknowns,  then  for  every  primal  Simplex 
basis  5  of  F  the  solutions  of  Bxi  =  fej  and  x^B  —  bJ  are  also  integer. 

Property  (1)  allows  very  efficient  execution  of  Simplex  iterations,  property  (2)  mo¬ 
tivates  the  use  of  special  data  structures  which  allow  efficient  storage  and  update  of 
Simplex  bases  and  property  (3)  allows  all  computatiors  to  be  performed  in  integer 
arithmetic  (assuming  ir  and  6  are  integer). 

Several  contributions  have  been  made  to  solving  variations  of  (P.N'SC)  when 
the  F-type  constraints  are  not  vacuous.  The  e  approaches  have  their  roots  in  work 
by  Haul  [1965]  and  Bennett  [1966].  Chen  and  Saigal  considered  a  version  of  (P.N'SC) 
with  all  equality  constraints.  Hartman  and  Lasdon  [1972]  considered  a  specialization 
of  (P.N’SC)  to  the  multicommodity  capacitated  transshipment  problem  (MCTP),  in 
which  the  F-type  constraints  are  generalized  upper  bound  (GUB)  constraints  and 
the  F-type  constraints  decouple  into  independent  pure  network  subproblems.  In 
each  of  these  treatments,  the  resulting  algorithm  is  a  specialization  of  the  primal 
Simplex  .Method.  McBride  [1972]  considered  (PN'SC)  in  the  context  of  the  mutual 
primal-dual  method,  with  all  constraints  assumed  to  be  equalities.  Each  of  these 
approaches  allows  a  basis  representation  which  may  be  factored  into  two  compo¬ 
nents:  an  explicit  part  of  dimension  rn  bv  rn  and  a  factored  part  of  dimension  p  by 
P- 

Our  inequality  formulation  of  (P.N'SC)  allows  a  dynamic  basis  representation 
where,  just  as  in  the  Cl'B  specialization,  both  the  dimension  of  the  factored  part 
and  the  dimension  (T  the  exjdicit  part  may  vary  from  one  iteration  to  the  next. 

Ill 


B.  THE  FACTORED  TABLEAU 


Recall  from  Chapter  3  the  primal  row  basis  at  an  arbitrary  point  in  the  solution 
process  is; 


(/)  1 

f  Fii 

F,2 

F,3  \ 

(22) 

F„ 

F,2 

Fi3 

(7.1) 

(222)  ' 

V  0 

0 

-7  / 

The  corresponding  nonbasic  constraint  rows  are  then: 


(0  /  -/  0  n  \ 

_  (”)  ^ 21  ^22  F 23 

(nb)  0-/0 
(21’)  \  E21  E22  /■-23  / 

and  the  principal  part  of  the  factored  tal)leau  is: 


(7.2) 


(..)  f::  C' CV  -  f:,  F,p  T:,  -  f:,  fC  F.,  +  fn  F,?  (  f  e  -  Fn  CV  F,0 

-F;.F,-,'F,;.ir.'F,,F:V  -F:;.i,‘,'(F,3  -  F,*,' F.j) 

(ml  .Cl'  -Ct'FiiF,','  .-iiV'Fn  -  FmF|,'F,3) 

'H')  -Fjj.Ci'  *  FiiFit  Fu.ir,'  F.,  bVf„F,-,'  -  F:,F,-,'  F;3  -  F.p F,,  Ci'' F,i  -  F,,F,VF;3i 

.  -FjiFFi'Fij.'n  FiiF,','  -F:iF,",'F|3  -  Fij.-iii'l Fn  -  FiiF,,  Fn) 


F.yFjj-  F33.hi''Fii  FiiFu  Ff 


~  C.2|/  U 


(7.3) 

The  F-type  constraints  represent  the  node-arc  incidence  matrix  of  a  graph.  For 
notational  convenience,  let  us  define  the  corresponding  directed  graph  explicitly  as 
Q  —  {AC  A]  where  C and  A  are  as  jweviously  defined.  Note  that  while  we  ch'^'^'^c 
to  interfiret  all  columns  of  ( I’.NSC)  as  arcs  in  C.  any  number  of  these  columns  may  be 
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null  with  respect  to  the  F-type  constraints  and  thus  may  be  interpreted  graphically 
as  null  arcs. 

Since  the  rows  of  [Fn  Fu  F13]  in  the  primal  row  basis  B  form  a  subset  of  the 
rows  of  F,  [Fii  Fi2  F13]  may  itself  be  interpreted  as  the  node-arc  incidence  matrix 
of  a  graph  Qb  =  {AsMb}-  is  a  subset  of  A'*,  but  in  general  Ab  ?!  A.  To 
see  this,  define  Ad  =  A^  —  As  and  consider  an  arc  a’  =  (n,. n^)  where  n,  6  Ab 
and  Tij  G  A^o-  Then  a’’  “spans”  the  partitioning  of  A"  into  A'b  and  Ad-  While  a’ 
is  a  doubleton  arc  with  respect  to  Q,  it  is  a  singleton  arc  with  respect  to  Qb-  We 
call  such  an  arc  a  “dynamic  singleton”,  and  we  will  discover  that  such  arcs  require 
special  handling  in  our  implementation. 

Since  a  nonsingular  Fu  must  exist,  it  is  well  known  that  the  columns  of  Fn 
represent  a  rooted  spanning  tree  defined  over  the  nodes  of  Ag,  and  that  Fn  may  be 
placed  into  upper  triangular  form  by  row  and  column  permutations.  If  we  choose  to 
represent  Fn  rather  than  Fn*  in  our  implementation,  we  are  interested  in  performing 
two  fundamental  operations: 

1.  Solving  linear  systems  of  the  type 


Fn^i  =  Ih 


and 


„2 


Fu 


where  cj  and  Z2  are  unknown  and  61  and  62  are  rationals  (not  necessarily 
integers),  and 
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2.  l-’lacing  Fn  in  upper  triangular  form,  where  results  from  a  column  ex¬ 
change,  a  row  exchange,  a  column  and  row  addition  or  a  column  and  row 
deletion  performed  on  Fn- 

The  literature  on  (PXE)  demonstrates  that  with  a  proper  choice  of  data  structures, 
operation  (1)  may  be  performed  very  efficiently  when  6i  and  62  are  integer,  and  that 
(2)  may  be  done  efficiently  when  the  claiss  of  updates  is  limited  to  column  exchanges. 
We  will  extend  these  existing  approaches  to  deal  with  (1)  and  (2)  in  the  (PNSC) 
setting. 

C.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COIVLMON  SUB¬ 
EXPRESSIONS  IN  COLUMN  GENERATION 

To  examine  the  actions  required  by  Generate_Column(LC),  suppose  we 
want  to  generate  column  LC  of  the  principal  part  of  the  tableau  (T..3)  and  place 
the  results  in  the  vector  =  (-1 , ^3- -4)^  which  is  partitioned  to  conform  to 
(7.3).  Rewriting  (7.3)  in  a  more  convenient  manner,  we  have: 


-r-.:  v,  -  n. 


:Tfi  1  -£7;  -£7; 


-•  11  *  •  u 

--•ij';'  £i:£;V 

-£7;  -FrrF,:  ^ 


~  ^5;  -^13  ~  j:  -  i 

i 

-O:  qv  V-u  -  f::  | 

F-rr,  ~FFF:  j 

■  r,,  -  £,:  r .*';3  ■ 

'a.FF,-  E„riF:,  ' 

-f.,  .-iy  ' 


The  calculations  may  then  proceed  as  follows: 


1.  Calculate  region  (in): 

(a)  If  column  LC  is  in  (j),  then: 

(b)  If  column  LC  is  in  (jj),  then: 

^3  =  -a;! 

(c)  If  column  LC  is  in  (jjj),  then; 

3,  =  .4.-.'  -  £.l^Tl'(f.3)‘'^] 

notice  that  and  may  be  computed  by  backpath  traversal 

on  Fn  (^-gM  Bradley,  Brown  and  Graves  [1977],  p.  II). 

2.  Calculate  region  (i): 

(a)  If  column  LC  is  in  (j),  then: 

Fll^l  =  Fi2^3 

(b)  If  column  LC  is  in  {jj),  then: 

Fn-^i  —  Fi2^3  + 

(c)  If  column  LC  is  in  {jjj),  then: 


F„m  =F,2r3  +  (Fi3)"^ 
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3.  Calculate  region  {ii): 


(a)  If  column  LC  is  in  (j),  then: 

Z2  =  — F22Z3  —  ^21^1 

(b)  If  column  LC  is  in  (jj),  then: 

Z2  =  —F22Z3  —  F2\Zi 

(c)  If  column  LC  is  in  (jjj).  then; 

Z2  =  —^'22-3  ~  +  (-^23)^^ 


4.  Calculate  region  (zr); 

(a)  If  column  LC  is  in  (j),  then: 

^•4  =  —F22Z2  —  E21Z1 

(b)  If  column  LC  is  in  {jj),  then: 

24  =  —E22Z3  —  E21Z1 

(c)  If  column  LC  is  in  {jjj),  then: 

24  =  —E22Z3  —  E21Z1  +  (£'23)^*" 

Then  column  LC  of  (7.4)  is  available  as  z^  =  {■21,22, 23, 24)^.  Note  that  floating 
point  multiplications  and/or  divisions  are  necessary  only  for  computation  of  regions 
(Hi)  and  (iv),  and  that  regions  (i),  (n)  and  (iv)  are  linear  combinations  of  region 
(lii)  and  one  another. 


D.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 


EXPRESSION  IN  ROW  GENERATION 

Now  we  examine  the  actions  of  Gencrate_Row{LR).  Suppose  we  desire  to 
generate  row  LR  and  place  the  results  of  the  computation  in  the  vector  z  =  (£5,  ig,  ir) 
which  is  partitioned  to  conform  to  (7.3).  Rewriting  (7.3)  in  a  convenient  form: 


(•) 


(n.) 

(tv) 


-  (i-f'i' Ai-'n':  -  rj  F-',‘ 

■  \[(  ^JI^n'r ij  -  £11  -  fji)  £1"' 


-  i4r,'i  €„Fr-' 
l£„fr,'r„  -  -  r,,)  4;,'; 


-Fv.'f.iAt; 
(£ii£iV£ It  “  £::'4~i 

4r! 


[“£11  '  £is  —  l~  1:  £'t''^ii  £n  —  £n  i  £13 

!_(£*:£ 1 1  ’  £1 1  ~  £:: )  4 1 1' ;  £i  s 

■t'  j-  {|(£i,£r,'£.:  -  £n)4rd  £u  ~  £;, }  f,*';  £.j  *  £n 
[4r.:i£„-{-j4r,''f„  £,-'}£.! 
jl  £:i£n'£ i)  “  £?:,(  4,/  i  £ij 
I-  (i(£„f7,'£i:  -  4r,'i  £„  ♦  £:, )  £,1'’  £|J  -  £0 


(7.5) 

Then  we  may  proceed  as  follows: 

1.  Calculate  region  (j): 

(a)  If  row  LR  is  in  (?),  then: 

h  =  -{Fu^)LRFn{Arl) 

(b)  If  row  LR  is  in  (if),  then: 

h  =  iF,,F~'F,,-F22)LRAn' 
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(c)  If  row  LR  is  in  (n?)>  then: 


is  =  iA^i)LR 

(d)  If  row  LR  is  in  (lu),  then: 

‘5  =  {^21^ ~  £22)lrAx^ 

2.  Calculate  region  [jj): 

(a)  If  row  LR  is  in  (f),  then: 

=  "-s^^n  ~  ^LR 

(b)  If  row  LR  is  in  (ii),  then: 

~6^n  =  ~^hE\\  —  {E2\)lr 

(c)  If  row  LR  is  in  [Hi),  then; 

hEn  =  -ism’ll 

(d)  If  row  LR  is  in  [iv).  then; 

^^E\\  =  —Z'iEw  —  {E2\)lr 

3.  Calculate  region  [jjj)- 

(a)  If  row  LR  is  in  (i),  then; 

~7  =  ^iEi3  +  ~6Ei3 

(b)  If  row  LR  is  in  [ii).  then: 

-T  =  ^s^^is  +  -gFia  +  {E23)lr 
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(c)  If  row  LR  is  in  (ni),  then: 


Z7  —  ZsEiJ  +  Z(iF\3 

(d)  If  row  LR  is  in  (zn),  then: 

h  =  +  Z^Fii  +  {Exi)LR 

The  row  LR  of  (7.5)  is  available  as  i  =  (25,26)57)-  Note  that  floating  point 
multiplications  and/or  divisions  are  avoided  for  terms  involving  F22F13  and  F23, 
and  that  regions  {jj)  and  (jjj)  are  linear  combinations  of  region  (/). 

E.  DATA  STRUCTURES 

Several  suites  of  data  structures  have  been  proposed  for  algorithms  designed 
to  solve  (PNE)  (Johnson  [1966],  Srinivasan  and  Thompson  [1973],  Glover,  Karney, 
Klingman  and  Napier  [1974],  Bradley,  Brown  and  Graves  [1977]).  Our  implementa¬ 
tion  is  based  on  the  last. 

The  purpose  of  the  data  structures  is  to  define  a  graph  that  represents  Fn. 
Such  a  graph  forms  a  rooted  spanning  tree  defined  over  the  nodes  of  Ag,  which 
we  denote  as  Qfh-  Note  that  because  Fh  is  nonsignular,  Fn  must  contain  at  least 
one  singleton  column.  Since  Fn  is  maintained  in  triangulated  form,  we  associate 
with  each  row  IR  the  column  JC  in  which  the  element  appearing  on  the  diagonal  of 
IR  occurs.  This  association  is  analogous  to  the  GUB  row/key  variable  association 
defined  for  the  GUB  algorithm.  We  represent  this  cissociation  with  the  array  KEY(), 
of  type  INTEGER  and  dimensioned  from  1  to  (rr?  +  n),  which  records  for  each  row 
of  Fii  the  extrinsic  index  of  the  corresponding  key  variable,  and  for  each  column  of 
Fii  the  extrinsic  index  of  the  corresponding  basic  factored  row.  KEY()  is  undefined 
for  other  row  and  column  indices. 
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We  require  knowledge  of  the  row  ordering  of  Fw  to  define  a  triangulation. 
P0()  is  an  INTEGER  array,  dimensioned  from  1  to  p,  which  for  each  row  IR  of  Fu 
records  the  extrinsic  index  of  the  row  that  follows  IR  in  the  current  triangulation. 
This  ordering  is  referred  to  as  preorder,  and  the  successor  of  IR  is  called  its  preorder 
successor.  P0()  is  undefined  for  factored  rows  not  currently  in  Fu-  Note  that  the 
triangulated  column  ordering  of  Fu  may  be  deduced  from  P0()  and  KEY(). 

P0()  and  KEY()  capture  the  triangulated  row  and  column  ordering  informa¬ 
tion  for  a  representation  of  Fu  but  provide  no  means  for  interpreting  this  repre¬ 
sentation  as  a  rooted  spanning  tree.  The  predecessor  function  p()  is  a  well  known 
method  for  representing  trees  and  thus  rooted  spanning  trees.  To  define  p{ik).  we 
associate  with  each  node  4  of  Qfh  the  row  index  of  the  offdiagonal  element  in  the 
column  of  Fu,  when  the  rows  are  ordered  to  correspond  to  a  triangulation  of 
Fu-  Graphically,  p(i)  may  be  iterated  recursively  to  trace  the  backpath  from  node 
i  to  the  distinguished  root  node  of  Of,,- 

We  represent  p()  by  the  INTEGER  array  P(),  dimensioned  from  1  to  p,  which 
for  each  row  IR  in  Fu  records  the  extrinsic  index  of  the  predecessor  of  IR.  It  is 
convenient  to  keep  a  record  of  the  sign  of  the  diagonal  element  of  each  row  of  Fu. 
We  use  the  sign  bit  of  P()  to  do  so.  Assume  node  IRP  is  the  predecessor  of  node 
IR  in  the  current  representation  of  Fu.  If  the  diagonal  element  in  row  IR  is  a  plus 
one,  then  P{IR)  =  IRP,  while  if  the  diagonal  element  in  row  IR  is  a  minus  one, 
P{IR)  =  —IRP-  In  the  graph  Of,,-,  this  may  be  interpreted  as  follows:  if  the  actual 
orientation  of  arc  =  KEY(IR)  is  the  same  a.s  that  recorded  in  Of,,-,  which  is 
{I R,  I RP).  then  P{IR)  >  0).  If  the  actual  orientation  of  arc  is  opposite  that 
recorded  in  Of,,^  then  P{IR)  <  0. 

The  final  data  structure,  D{),  is  an  INTEGER  array  dimensioned  from  1  to  p. 
For  each  node  IR  of  Of,,  (row  of  F,,)  D(IR)  records  the  depth  of  node  IR,  where 


depth  is  defined  as  the  number  of  nodes  encountered  on  the  backpath  between 
IR  and  the  root  node.  D()  allows  updates  to  be  performed  on  the  data  structures 
representing  Qfu  ^  single  pass  through  the  data  structures,  and  allows  us  to  avoid 
unnecessary  operations  in  the  solution  of  linear  systems  to  be  discussed  shortly. 

F.  SOLVING  LINEAR  SYSTEMS 

It  is  clear  from  the  discussion  of  Generate_Row(LR)  and  Generate.Col- 
umn(LC)  that  the  solution  of  the  linear  systems, 


and 

where  and  :2  ^re  vectors  of  unknowns  and  bi  and  62  are  vectors  of  rationals,  are 
crucial  steps.  The  practical  value  of  our  implementation  depends  to  a  large  extent 
on  the  efficiency  with  which  these  systems  may  be  solved,  and  thus  we  will  consider 
this  problem  in  some  detail. 

First,  we  examine  the  data  structures  used  to  represent  the  terms.  Three  data 
structures  are  used  to  represent  61  and  b^,  the  right-hand  side  terms.  \VORK()  is  a 
DOUBLE  PRECISION  array,  dimensioned  from  1  to  p,  which  is  used  to  hold  the 
real  values  of  61  and  b^  (we  sequence  operations  so  that  at  any  given  moment  we  are 
interested  in  either  61  or  62  but  not  both,  so  that  the  same  array  may  be  used  for 
both).  An  INTEGER  array  \VORKMK(),  also  dimensioned  from  1  to  p,  is  used  as 
a  nonzero  ma^k  for  VVORK().  The  convention  is: 


WORKMK(K)  =  0  if  WORK(K)  =  0.0 
WORKMK(K)  =  1  IF  WORK(K)  ^  0.0 


Finally,  the  array  WORKNZ(),  of  type  INTEGER  and  also  dimensioned  from  1 
to  p,  records  the  extrinsic  column  (in  5i)  or  row  (in  62)  index  of  each  nonzero  in 
WORK().  The  use  of  WORKMK()  and  WORKNZ()  allows  the  efficient  manage¬ 
ment  of  nonzeros  in  bi  and  b^. 

The  solution  vectors  Zi  and  Z2  are  represented  by  three  analogous  arrays. 
ROWCOL()  (for  rowcolumn())  is  of  type  DOUBLE  PRECISION  and  is  dimensioned 
from  1  to  (m  +  n).  It  is  used  to  store  the  representation  of  a  row  and  column  of  the 
principal  part  of  the  tableau,  and  thus  portions  of  it  are  used  to  store  the  solutions 
of  Equation  (7.1).  Associated  with  ROWCOL()  is  ROWCOLMK(),  an  INTEGER 
array  of  conformable  dimension,  used  as  a  nonzero  mask  for  ROWCOL(),  and  the 
INTEGER  array  ROWCOLNZ(),  also  of  conformable  dimension,  used  to  record  the 
indices  of  nonzeros  in  RO\VCOL(). 

In  contrast  to  the  usual  situation  in  Simplex-based  approaches  to  linear  pro¬ 
gramming  in  which  each  successive  right-hand  side  of  the  problem  may  be  computed 
by  means  of  a  simple  update  to  the  previous  right-hand  side,  in  this  setting  such 
is  not  the  case.  Thus,  we  are  not  able  to  derive  the  solution  to  the  current  form 
of  Equation  (7.1)  by  simply  updating  the  previous  solution.  Instead,  we  must  solve 
each  system  from  scratch. 

Each  right-hand  side  is  itself  the  result  of  a  sequence  of  preliminary  computa¬ 
tions.  As  each  new  nonzero  element  is  generated  by  these  computations,  its  value 
is  placed  in  \VORK(),  \VORKMK()  is  updated  and  the  index  of  VVORK()  in  which 


the  nonzero  appears  in  placed  in  the  next  available  position  in  \VORKNZ().  For  ex¬ 
ample,  suppose  the  first  nonzero  generated  for  VVORK()  is  the  value  6.5  in  position 
38.  Then  the  arrays  appear  as  follows: 

WORK(38)  =  6.5 

\VORKMK(38)  =  1.0 

VV0RKNZ{1)  =  38.0 

A  counter  for  the  number  of  nonzeros  in  WORK()  is  maintained  (an  INTEGER 
variable  called  NNZVV),  so  at  the  conclusion  of  the  preliminary  computations  we 
may  iterate  through  the  nonzeros  on  the  right-hand  side  in  the  order  in  which  they 
were  generated  by: 


\VORK(\VORKNZ(K)),  K  =  1,  ...,  NNZVV 

Since  our  representation  of  Fn  is  in  upper  triangular  form,  the  obvious  ap¬ 
proach  for  solving 


—  ^2  (~-6) 

is  by  back  substitution.  .All  nonzero  elements  in  Fn  are  either  plus  ones  or  minus 
ones,  so  only  assignments,  additions  and  subtractions  are  necessary.  Assuming  that 
the  current  dimension  of  Fn  is  k  by  k,  we  proceed  by  considering  each  row  IR  in 
turn,  IR  =  k,....\.  We  assign  the  value  in  \VORK(IR)  to  RO\VCOL(KEV(IR)) 
and  update  by  an  addition  the  value  in  W'ORK()  in  the  row  corresponding  to  the 
predecessor  of  row  IR.  (i.e.,  \VORK(P(IR))). 

This  approach  requires  knowledge  of  the  ordering  of  the  rows  of  Fn  in  the 
order  of  last  to  first  in  triangulated  form,  which  is  precisely  the  reverse  of  our  data 
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Figure  7.1:  A  Triangulated  Fn 

structure  PO().  To  support  this  approach  in  our  implementation,  we  could  define 
and  maintain  an  additional  data  structure  E0(),  recording  the  “endorder”  of  Fn.  If 
we  choose  not  to  include  an  additional  data  structure,  we  could  instead  invert  P0() 
in  situ  whenever  this  reverse  ordering  is  required,  inverting  it  again  when  P0()  is 
next  required.  With  either  approach,  it  is  convenient  to  make  P0()  and  E0().  if 
included,  circular  lists  and  to  add  a  distinguished  “artificial  node”,  say  IRA,  which 
is  then  used  to  locate  both  the  first  node  (row)  in  preorder  and  the  first  node  in 
endorder.  To  illustrate.  Figure  7.1  displays  a  triangulated  form  of  Fn  with  row  and 
column  labels. 

Figure  7.2  presents  the  corresponding  graph  By  assuming  the  existence 
of  an  artificial  node  IR.A.  we  in  effect  change  our  paradigm  of  Qfu  '^hat  shown  in 
Figure  7. .3. 

We  may  interpret  the  backsolve  method  as  a  labeling  procedure  on  In 

this  context,  each  nonzero  in  the  right-hand  side  is  interpreted  a.s  a  supply  or  demand 
at  the  corresponding  node  in  Gf^^-  The  problem  of  solving  Equation  (7.6)  is  then 
interpreted  as  one  of  finding  a  set  of  feasible  flows  defined  on  the  arcs  of  Qf^- 

1 21 
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Figure  7.2:  A  Basis  Graph  (Jf  ,, 

The  backsolve  approach  to  solving  Kfpiation  (7.1)  lias  tlie  advantage  of  sim¬ 
plicity,  and  when  viewed  a.s  a  labeling  algorithm  has  the  intuitive  appeal  of  “visiting” 
each  node  only  once.  The  disadvantages  are  the  need  for  either  additional  storage 
(E0())  or  additional  computation  (two  inv<-rsions  of  P0()). 

The  right-hand  sides  of  these  probh'ms  are  invariably  \ery  sjjarse;  the  density 
is  typically  0..3'T  or  less.  We  are  thus  motivated  to  exjilore  alternate  approaches  to 
solving  Equation  (7.6).  I  he  solution  approach  just  outlined  recpiircs  the  “visiting"  of 
every  node  in  Qf-., .  Although  W’OHKMK()  allows  us  to  perform  a  simple  l.NTEGER 
comparison  to  determine  if  the  current  node  has  nonzero  supply  or  demand,  we  have 
strong  empirical  evidence  that  suggests  (hat  other  approaches  are  more  efficient.  It 
is  important  to  kcx'])  in  mind  that  tliis  problem  is  very  deeply  nested  within  the 
algorithm  st  nu  t  tire,  and  must  be  stdved  tfiis  of  thousands  of  times  for  a  typical 


Figure  7.3:  An  Alternate  Graph  Paradigm 

linear  programining  problem.  Small  changes  in  efficiency  here  have  large  effects  on 
the  overall  algorithm  efficiency. 

Since  tlie  right-hand  sides  are  sparse,  an  alternative  is  to  consider  solving  Equa¬ 
tion  (7.6)  iteratively  as  a  sequence  of  subproblems,  each  consisting  of  a  right-hand 
side  containing  a  single  nom'cro  element  and  a  current  partial  solution.  Iterating 
through  the  nonzeros,  we  introduce  the  supply  or  demand  at  the  corresponding  node 
in  ,  and  then  adjust  the  flows  of  all  arcs  appearing  on  the  backpath  of  that  node 
as  necessary. 

Analyzing  the  worst-case  performance  of  these  two  approaches,  we  see  that  the 
first  approach  is  linear  in  the  number  of  nodes  in  (rows  in  Fn),  while  the  second 
approach  is  quadratic  in  the  number  of  nonzeros  in  the  right-hand  side.  For  the 
densities  typically  encountered  in  the  problems  we  have  studied,  the  crossover  point 
at  which  the  performance  of  the  linear  algorithm  overtakes  that  of  the  quadratic  is 


126 


in  the  range  of  400  to  1000  rows  in  Fu-  However,  in  the  testing  we  have  performed, 
the  worst-ca^e  analysis  is  pessimistic,  and  in  practice  the  performance  of  the  second 
approach  always  dominates  that  of  the  first. 

The  second  linear  system  of  interest  is 

zjFn-^bJ  .  (7.7) 

We  interpret  the  problem  as  one  of  being  given  flows  on  the  arcs  of  and  being 
a^ked  to  find  the  necessary  supplies  and  demands  at  the  nodes.  Analogies  to  each  of 
the  approaches  to  solving  Equation  (7.6)  may  be  found  for  solving  Equation  (7.7), 
and  our  empirical  evidence  strongly  suggests  that  an  approach  analogous  to  the 
second  method  for  solving  Equation  (7.6)  is  preferable  for  solving  Equation  (7.7). 

Having  made  the  design  decisions  for  solving  Equation  (7.6)  and  Equation 
(7.7),  it  is  worthwhile  to  review  our  paradigm  for  Qfm-  Our  preferred  solution 
techniques  do  not  require  a  partial  ordering  of  all  rows  of  Fu.  Rather,  we  require 
partial  ordering  information  only  among  each  set  of  coupled  rows,  or,  graphically, 
among  each  disjoint  component  of  Qfh-  Rather  than  introducing  an  artificial  node, 
we  may  treat  each  disjoint  component  as  an  independent  entity.  It  is  then  conve¬ 
nient  to  treat  singleton  arcs  (either  dynamic  or  static  singletons)  as  self-loops,  and 
construct  the  predecessor  function  so  that  it  forms  a  ring  within  each  component. 
Our  graph  paradigm  of  Figure  7.1  then  becomes  as  shown  in  Figure  7.4.  The  data 
structure  representation  of  Figure  7.4  is  then; 
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Figure  7.4:  The  Implementation  Paradigm  for 

G.  FACTORED  KERNEL  UPDATE 

The  update  of  a  rooted  spanning  tree  representation  of  a  pure  network  tri¬ 
angulated  basis  wliicli  results  from  a  column  exchange  is  well  understood  and  has 
been  treated  thoroughly  in  the  literature.  Bradley,  Brown  and  Graves  [1977]  give 
an  excellent  account  of  such  an  update  in  an  algorithmic  setting  very  similar  to  the 
one  here.  Hence,  we  will  not  repeat  the  details.  In  this  setting,  however,  a  column 
exchange  is  but  one  of  four  possible  updates  required.  We  must  also  deal  with  row- 
exchanges,  row  and  column  additions  and  row  and  column  deletions.  Our  general 
approach  will  be  to  reduce  these  three  cases  to  the  column  exchange  case  through 
a  sequence  of  operations  which  may  be  thought  of  as  pre-  and  post-processing.  The 
following  three  cases,  along  with  the  column  exchange  case,  comprise  the  actions 
required  for  Update JFactored_Kernel. 

The  first  case  we  consider  is  the  row  and  column  deletion  case,  in  which  the 
dimension  of  Tu  decreases.  It  is  convenient  lo  limit  row/column  deletions  to  row/key 
variable  pairs,  d  his  is  because  removing  a  row/key  variable  pair  always  preserves 
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nonsingularity  of  the  remaining  factored  kernel,  and  because,  with  our  choice  of  data 
structures,  given  one  member  of  the  pair,  it  is  ea^y  to  identify  the  other  member 
through  the  use  of  the  KEY()  array. 

With  our  paradigm  for  Qpu  >  hey  variable  for  row  IR  corresponds  to  the  arc 
which  connects  node  IR  to  its  predecessor.  The  removal  of  IR  and  KEY(IR)  creates 
a  new  disjoint  component  within  consisting  of  all  nodes  and  the  associated  key 
vajiables  on  whose  backpaths  IR  and  KEY(IR)  appear.  If  IR  is  a  leaf  (a  node  which 
is  incident  to  a  single  arc),  then  the  new  disjoint  component  is  null.  Within  the  new 
disjoint  component,  D()  changes  for  all  nodes,  but  P()  and  P0()  change  only  for  the 
first  node  in  preorder,  say  IRFIRST,  and  the  last  node  in  preorder,  say  IRLAST. 
For  the  first  node  in  preorder,  its  predecessor  becomes  itself  (with  the  sign  bit  used 
to  indicate  arc  orientation),  forming  a  self-loop.  The  update  is  thus: 

P  (IRFIRST)  =  IRFIRST 

Its  depth  is  updated  to  be  zero,  and  the  magnitude  of  this  depth  change  is  recorded 
for  future  use.  For  every  other  node  in  the  new  disjoint  component,  D{)  is  reduced  by 
the  magnitude  of  the  IRFIRST  depth  change.  Finally,  PO()  of  IRL.AST  is  changed 
to  the  first  node  in  preorder,  restoring  the  local  ring  structure  of  P0().  The  update 
is: 


PO  (IRLAST)  =  IRFIRST 

These  updates  can  be  accomplished  in  one  pass  through  the  data  struciures  of  the 
new  disjoint  component.  The  structure  of  Qfu  is  restored  by  changing  P0()  of  the 
predecessor  of  IR  to  the  index  of  the  node  which  was  the  successor  of  IRL.AST  prior 
to  the  update. 
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The  second  update  case  is  a  row  and  column  addition.  Suppose  IRADD  is  the 
row  (node)  to  be  added  to  Fn  (^Fh)-  wish  to  treat  this  update  as  a  column 
exchange.  To  do  so,  we  first  incorporate  IRADD  into  the  structure  and  associate 
with  it  an  imaginary  (“logical”)  arc  which  forms  a  self-loop.  We  then  perform  a 
column  exchange  with  the  new  column,  say  JCADD,  entering  Qfu  logical 

arc  leaving  Gfi,- 

The  difficulty  in  incorporating  IRADD  into  Qfh  lies  with  the  dynamic  single- 
ton  arcs.  Qfu  currently  may  have  many  singleton  arcs  which  appear  in  our  paradigm 
cis  self-loops.  If  any  of  these  dynamic  singletons  are  actually  incident  to  IRADD  in 
G,  incorporating  IRADD  into  Gpn  requires  changing  their  representation. 

To  illustrate  this  point,  assume  IRADD  has  the  node  label  “9”  and  we  wish 
to  incorporate  it  into  the  graph  shown  in  Figure  7.4.  Assume  also  that  arcs  103  and 
107  are  dynamic  singletons  which  are  actually  incident  to  node  9;  that  is,  arc  103 
=  (9,4)  and  arc  107  =  (9,8).  We  first  initialize  the  data  structures  for  node  9  as 
follows: 


D(9)  =  0 
F(9)  =  9 
PO{9)  =  9, 

which  has  the  effect  of  placing  node  9  at  depth  0  as  a  disjoint  component  consisting 
of  a  single  node  with  a  (logical)  self-loop,  as  shown  in  Figure  7.5. 

We  then  change  the  representation  of  arcs  103  and  107  from  dynamic  singletons  to 
doubletons.  To  do  this,  we  change  P(4)  from  -4  to  -9  and  P(8)  from  -8  to  -9.  We 
then  assert  a  partial  ordering  among  those  (formerly)  disjoint  components  which 
have  been  merged  by  the  change  in  the  representation  of  the  dynamic  singletons. 
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Figure  7.5:  Initialization  of  Node  9  Prior  to  Incorporation  into 

and  retriarigulate.  'I'his  is  acconiplislicd  hy  increasing  the  depth  of  each  node  in  tlie 
component  by  one,  changing  PO(7)  from  1  to  8,  and  PO(8)  from  8  to  9.  This  update 
requires  a  single  pass  through  the  data  structures  of  each  affected  component.  After 
incorporation,  appears  as  shown  in  figure  7.6.  We  may  now  complete  the 

update  by  performing  a  column  exchange. 

To  identify  the  dynamic  singletons  which  are  affected  by  such  an  update,  we 
access  the  extrinsic  problem  data  structures  for  row  (node)  IRA  DU.  For  each  column 
with  a  nonzero  in  row  IRADD,  we  test  whether  or  not  that  column  is  currently  key 
for  some  row  in  Fn-  If  so,  that  column  is  currently  represented  in  as  a  dynamic 
singleton. 

The  final  update  case  is  the  row  exchange  case.  Suppose  we  want  to  replace 
row  IROUT  with  row  IRIN.  Wc  treat  this  in  two  stages.  The  first  stage  is  a  row 
and  column  deletion  case,  in  which  the  node  IROUT  and  the  arc  KE'j  (IROUT)  are 
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removed  from  Qfh-  To  complete  the  update,  we  seek  a  column  (arc)  which  may 
be  added  to  along  with  row  (node)  IRIN.  We  require  an  arc  which  is  either 
a  static  singleton,  incident  to  node  IRIN,  or  a  doubleton,  incident  to  node  IRIN 
and  a  node  (other  than  IROUT)  which  is  currently  in  Qfu-  We  first  consider  arc 
KEY(IROUT).  If  is  satisfies  the  second  condition  (obviously,  it  cannot  satisfy  the 
first  condition),  we  designate  it  as  the  arc  to  be  added.  If  not,  we  search  among  the 
variables  (arcs)  in  region  {Hi)  of  the  tableau  for  such  an  arc.  We  know  one  must 
exist,  for  otherwise  B  is  singular.  We  select  the  first  such  arc  found  as  the  arc  to  be 
added.  We  then  perform  the  row  and  column  addition  case. 

Recall  that  the  action  Is_Factored_Kernel_Singular()  is  required  when  we 
have  identified  an  arc  (column)  for  removal  from  Qf^  a-nd  we  are  considering  a 
candidate  arc  to  replace  it.  We  want  to  determine  if  this  exchange  preserves  the 
nonsingularity  of  Fu-  The  fact  that  Qfu  is  a  rooted  spanning  tree  provides  the 
necessary  structure  to  support  a  simple  test  of  nonsingularity.  Assume  JCOUT  is 
the  (known)  arc  to  be  removed  from  Qfu  ^"d  JCIN  is  the  candidate  replacement. 
In  our  paradigm,  Fu  is  nonsingular  if  and  only  if  each  disjoint  component  of  Qf^ 
is  connected  and  contains  exactly  one  cycle,  which  must  be  a  self  loop  occurring  at 
the  component’s  root  (the  root  is  the  distinguished  node  in  the  component  whose 
depth  is  zero).  The  removal  of  JCOUT  creates  a  new  disjoint  component  which  has 
no  root  self  loop.  If  the  addition  of  JCIN  fails  to  correct  this,  the  proposed  exchange 
is  singular. 

The  specific  test  is  this:  identify  the  nodes  incident  to  JCIN  using  the  original 
problem  data  structures.  Note  that  there  may  be  zero,  one  or  two  such  nodes.  If 
none  of  these  nodes  are  currently  included  in  the  structure  of  Qfu{Fi\).,  then  the 
proposed  exchange  is  singular.  For  each  such  node  which  is  currently  included  in 
the  structure  of  traverse  the  backpath  of  that  node  by  recursively  iterating 


the  predecessor  array  P{)  until  either  the  root  node  of  that  component  is  found,  or 
until  the  “join”  is  located  (the  “join”  is  that  node  with  largest  depth  which  appears 
on  the  backpath  of  both  nodes  incident  to  JCIN.  The  join  exists  only  if  JCIX  is 
incident  to  two  nodes,  both  nodes  are  currently  included  in  the  structure  of  Qfui 
and  the  two  nodes  are  in  the  same  disjoint  component.).  If  JCIN  is  encountered 
during  this  backpath  traversal,  the  proposed  exchange  is  nonsingular.  Otherwise,  it 
is  singular. 

The  action  Find_Column_to_Remove(IOUT)  is  particularly  simple,  since, 
as  mentioned  previously,  we  always  remove  row/key  variable  pairs.  Thus,  if  IR  is 
the  node  to  be  removed,  then  KEY(IR)  is  the  corresponding  column  to  be  removed. 

Finally,  the  action  Find_Column_to_Add(IIN)  requires  searching  among 
the  variables  in  region  {Hi)  of  the  tableau.  For  each  candidate  arc,  we  invoke 
IsJFactored_KerneLSingular().  If  the  response  is  “FALSE”,  we  have  found  the 
column  to  be  added.  Otherwise,  we  continue  the  search. 


133 


VIII.  FACTORIZATION  OF  GENERALIZED 

NETWORK  ROWS 

A.  INTRODUCTION 

The  problem  of  interest  remains: 

min  :  wy 
y 

(GNSC)  s.t.  :  Fy  <  6;  F  phyn 
Ey  <  r;  E  m  by  n 
—  ly  <0;  —In  by  n 

where  F,  F,  —7,  w,  b  and  r  are  as  before.  We  now  require  that  each  column  of 
F  have  at  most  two  nonzero  elements,  w'hich  may  be  of  arbitrary  sign.  We  may 
associate  a  generalized  network  with  F  by  defining  a  node  corresponding  to  each 
row  i  of  F  and  an  arc  F-'  corresponding  to  each  column  j  of  F.  The  arcs  are  defined 
aLs; 

(?,  k)  if  F,j  ^  0,  Fjtj  ^  0  and  i  <  k  <  p 

F-'  =  -  (?, 0)  if  F,j  ^  0,  Fkj  =0  for  all  k  ^  i  <  p 

(0.0)  if  F,j  —  0  for  all  ^  <  p 

VVe  let  A  =  {F' . F"}  ,  A”  =  {1 , . . .  ,p}  and  define  the  graph  Q  =  {A^,  A}.  Arcs 

of  the  form  (i,0)  are  singleton  arcs,  sometimes  called  root  arcs.  Arcs  of  the  form 
(0,  0)  are  null. 

The  most  widely  studied  specialization  of  (GNSC)  is  obtained  when  the  F- 
type  constraints  are  equalities  and  the  F-type  constraints  are  vacuous,  resulting 
in: 
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min  :  wy 
y 

(GNE) 

s.t.  :  Fy  <  b;  F  p  hy  n 
Ey  <  r  E  mhy  n 
—  ly<0\  —I  n  by  n 

Dantzig  [1963]  provides  the  seminal  treatment  of  (GNE),  identifying  the  im¬ 
portant  structure  that  leads  directly  to  efficient  primal  Simplex-btised  algorithms 
for  its  solution.  Implementations  have  been  reported  by  Glover,  Klingman,  Hultz, 
Karney  and  Elam  [1972.  1973,  1977,  1978,  1979]  and  Brown  and  McBride  [1984]. 

(GNE)  may  be  viewed  as  a  generalization  of  (PNE),  and  the  cost  of  such 
generalization  is  a  weakening  of  the  properties  which  so  strongly  characterize  (PNE) 
and  which  lead  to  the  efficient  and  elegant  implementations.  In  (GNE),  F  is  not 
totally  unimodular,  and  thus  an  implementation  may  not  be  restricted  to  integer 
arithmetic.  It  is  not  possible  to  triangulate  every  primal  Simplex  basis  extracted 
from  F  by  row  and  column  permutations.  The  subgraph  generated  from  the  columns 
of  a  primal  Simplex  basis  no  longer  form  a  rooted  spanning  tree  defined  on  the  nodes 
of  Q.  Apparently,  much  has  been  given  up  in  the  generalization. 

In  fact,  significant  structure  re.nains  in  (GNE).  A  well  known  result,  due  to 
Dantzig  [1963].  shows  that  any  primal  Simplex  basis  extracted  from  F  can  be  put  in 
the  form  of  Figure  (8.1)  by  row  and  column  permutations.  Each  square  submatrix 
component  fi*'  is  either  upper  triangular  or  nearly  upper  triangular  with  only  one 
element  below  the  diagonal.  Notice  that  if  each  B*'  is  strictly  upper  triangular,  the 
structure  is  analogous  to  that  found  in  the  (PNE)  primal  Simplex  basis. 

Interpreting  the  structure  of  as  a  subgraph  of  Q,  we  find  that  when 
is  upper  triangular,  the  subgraph  may  be  viewed  as  a  disjoint  component  with  a 
single  selMoop  at  the  root  node,  exactly  as  in  the  (PNE)  case.  When  B*  is  not 
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Figure  8.1:  Near-Triangulated  Simplex  Basis  Corresponding  to  (GNE) 

upper  triangular,  the  subgraph  still  forms  a  disjoint  component  with  a  single  loop, 
but  that  loop  is  no  longer  a  self-loop.  For  example.  Figure  (8-2)  shows  a  nearly 
triangular  B*'  with  row  and  column  labels. 
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Figure  8.2;  A  Sample  Nearly-Triangulated  (GNE)  Simplex  Basis  Com¬ 
ponent 

The  resulting  subgraph  is  shown  in  Figure  (8.3).  Such  a  subgraph  is  commonly 
called  a  “one-tree”,  (an  apparent  oxymoron). 

This  structure  allows  the  extension  of  the  algorithm  and  data  structures  de¬ 
veloped  for  (PNE),  leading  to  efficient  implementations  to  solve  (GNE)  (e.g..  Brown 
and  McBride  [1984]). 

(GNSC)  has  received  less  attention  in  the  literature  than  (PNSC).  Hultz 
and  Klingman  [1976]  develop  a  primal  Simplex- based  approach  to  an  equality- 
constrained  formulation  of  (GNSC).  and  report  on  an  implementation  which  allows 
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Figure  8.3:  “One-tree”  Subgraph 

a  single  £-type  constraint  (Ilultz  and  Kltngrnan  (1978)).  McBride  [}98oj  develops 
an  algorithm  for  solving  a  generalization  of  (GNSC)  in  which  complicating  columns 
as  well  as  complicating  constraints  are  permitted,  and  reports  on  an  implementation 
of  that  algorithm. 

B.  THE  FACTORED  TABLEAU 

The  algebraic  development  of  the  primal  row  basi.s  B.  the  nonbasic  rows  D 
and  the  factored  tableau  is  exactly  as  shown  in  Chapter  7  and  is  not  repeated  here. 
Note  that  F]i,  the  factored  kernel,  may  now  be  placed  in  the  form  shown  in  Figure 
(8.1)  by  row  and  column  permutations.  The  corresponding  graph,  consists  of 
one  or  more  disjoint  components,  each  of  which  contains  exactly  one  loop.  The  loop 
may  be  either  a  self-loop,  as  in  the  case  of  the  disjoint  components  of  (PNSC),  or 
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it  may  be  a  "root  loop”  (a  loop  consisting  of  two  or  more  nodes,  all  of  which  are  at 
depth  zero  in  the  component),  as  shown  in  Figure  (8.3). 

C.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 
EXPRESSIONS  IN  COLUMxN  GENERATION 

The  sequencing  of  computations  for  column  generation  in  (GNSC)  is  exactly 
the  same  as  that  for  (PNSC). 

D.  SEQUENCING  COMPUTATIONS  TO  EXPLOIT  COMMON  SUB¬ 
EXPRESSIONS  IN  ROW  GENERATION 

The  sequencing  of  computations  for  row  generation  in  (G.NSC)  is  exactly  the 
same  as  that  for  (PNSC). 

E.  DATA  STRUCTURES 

Several  proposed  data  structures  for  algorithms  tailored  to  solve  (GNE)  have 
been  offered  (Glover,  Klingman  and  Stutz  [1973,  1974],  Elam,  Glover  and  Klingman 
[1979],  Brown  and  McBride  [1984]).  Our  implementation  is  based  on  the  last,  which 
extends  to  (GNE)  the  data  structures  developed  in  Bradley,  Brown  and  Graves 
[1977]  to  solve  (PNE). 

We  have  seen  the  similarity  between  the  structure  of  the  disjoint  components 
that  arise  in  the  graph  corresponding  to  Fn  in  (GNSC)  and  the  structure  of  the 
components  in  (PNSC).  To  develop  a  representation  of  these  components,  we  wish  to 
extend  the  data  structures  developed  for  (PNSC)  in  a  natural  way.  Knowledge  of  the 
row  ordering  of  Fn  is  again  maintained  in  an  INTEGER  array  PO{).  dimensioned 
fron;  1  to  p.  .As  before,  it  is  convenient  to  make  P0{)  a  circular  list  defined  on 
each  disjoint  component.  The  unique  correspondence  between  each  row  of  Fn  and 
a  distinguished  column  is  maintained  in  the  INTEGER  array  KE^'()-  dimensioned 


from  1  to  (m  +  n).  The  conventions  for  KE'N’O  are  exactly  as  in  (PNSC).  The  depth 
array,  D{),  is  defined  exactly  as  before.  However,  as  suggested  by  our  diagram 
in  Figure  (8.3),  we  adopt  the  convention  of  placing  all  nodes  appearing  on  a  loop 
at  depth  zero,  which  is  consistent  with  our  definition  of  "root  loop".  Finally,  a 
representation  of  the  predecessor  function  is  needed,  and  we  define  P{)  as  before. 
For  any  row  IR  in  the  factored  kernel,  P{J R)  >  0  indicates  that  the  diagonal  element 
in  the  near-triangulation  is  the  first  non-zero  factored  element  in  its  column. 

F.  SOLVING  LINEAR  SYSTEMS 

Just  as  in  (PNSC),  the  crucial  steps  in  the  generation  of  rows  and  columns  of 
the  principal  part  of  the  tableau  are  solving  the  systems 

=lFu=bJ  (8.1) 

and 


Fu=2  =  62  (8.2) 

where  Zi  and  C2  are  vectors  of  unknowns  and  and  b-2  are  vectors  of  rationals.  Our 
approach  for  solving  these  systems  closely  parallels  that  developed  for  (PNSC). 

ihe  data  structures  used  to  represent  c,.r2,6i  and  62  are  exactly  the  same 
as  in  the  (P.NSC)  implementation.  To  represent  bj  and  62,  we  use; 

\VORK() 

WORKMKO 


WORKNZO, 


and  for  and  Z2: 

ROWCOLO 

ROWCOLMKO 

ROWCOLNZ(). 

The  right-hand  side  sparsity-exploiting  approach  for  solving  (8.1)  and  (8.2) 
developed  in  Chapter  7  may  still  be  used,  but  two  complications  arise  in  the  (GNSC) 
setting.  First,  the  nonzero  elements  of  Fn  are  no  longer  restricted  to  be  plus  and 
minus  ones,  and  may  assume  arbitrary  values.  Thus,  we  are  not  able  to  restrict  all 
arithmetic  operations  to  additions  and  subtractions. 

Interpreting  (8.1)  and  (8.2)  as  problems  of  finding  balancing  supply  and  de¬ 
mands  (8.1)  or  feasible  flows  (8.2),  we  see  that  in  (GNSC)  we  must  allow  for  gains 
and  losses  in  the  flows  over  the  arcs  of  Because  our  formulation  allows  two  arbi¬ 
trary  nonzeros  in  each  column  of  F  (rather  than  one  arbitrary  nonzero  and  one  plus 
one,  for  example),  backward  and  forward  substitution  schemes  require  both  multi¬ 
plications  and  divisions.  W’e  may  eliminate  the  need  for  divisions  by  pre-computing 
both  ratios  of  the  nonzero  elements  in  each  arc.  That  is.  if  a  column  has  nonzero 
elements  a  and  b  in  rows  of  F,  we  pre-compute  the  ratios  |  and  We  then  substi¬ 
tute  the  pre-computed  value  for  the  division  whenever  it  occurs.  Storing  a  pair  of 
DOUBLE  PRECISION  real  numbers  for  every  column  of  (GNSC)  is  wasteful,  and 
instead  we  define  a  single  pointer  for  each  column  which  points  to  the  location  of 
the  first  of  two  DOEBLE  PRECISION  real  values.  The  first  is  the  value  |  and  the 
second  is  E  (Note  that  for  singleton  arcs,  the  values  a  and  ^  are  stored  instead.) 
We  compute  and  store  only  the  unique  ratio  pairs  for  a  given  problem  instance,  and 
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the  number  of  such  unique  pairs  is  tisually  quite  small:  typically,  a  problem  with 
10,000  to  20,000  columns  has  only  30  to  40  unique  ratio  pairs. 

The  second  difficulty  is  that  the  disjoint  components  of  Fn  may  be  nearly  tri¬ 
angular  rather  than  triangular,  and  thus  backward  and  forward  substitution  cannot 
be  applied  directly.  However,  Dantzig  [1963]  shows  that  a  variation  of  backward 
and  forward  substitution  may  be  used  to  soK'e  (8.1)  and  (8.2).  This  method  solves 
the  triangulated  part  of  a  disjoint  component  exactly  as  is  done  in  backward  and 
forward  substitution.  When  a  root  loop  is  encountered,  the  method  requires  two 
calculations  for  each  row  or  column  of  Fn  corresponding  to  a  node  or  arc  on  the 
loop.  The  second  calculation  involves  a  term  Dantzig  [1963]  called  the  “loop  factor"’, 
which  is  common  value  for  every  node  (or  arc)  on  the  loop.  Our  implementation 
uses  this  modified  forward  and  backward  solution  technique  for  solving  (8.1)  and 
(8.2),  and  we  store  unique  values  of  loop  factors  in  a  DOUBLE  PRECISION  array 
called  GNROOTO. 

G.  FACTORED  KERNEL  UPDATE 

We  first  consider  the  actions  required  by  UpdateJFactored_Kernel. 

Brown  and  McBride  [1984]  give  an  excellent  treat:. lent  of  the  update  required 
for  a  column  exchange  within  a  generalized  network  basis  in  an  algorithmic  setting 
very  similar  to  the  one  here,  and  thus  we  do  not  repeat  that  discussion.  As  in 
the  case  of  (PN'SC).  however,  three  additional  clcisses  of  updates  may  occur  in  this 
setting:  row  exchanges,  column  and  row  additions  and  column  and  row  deletions. 
We  again  treat  these  three  ca.ses  by  reducing  them  to  the  column  exchange  case 
through  sequences  of  pre-  and  post-processing  operations. 

We  again  limit  row  and  column  deletions  to  row/key  variable  pairs.  When 
row  IR  and  its  a.ssociatf^d  key  variable  KE^’(IR)  are  removed  from  Fn,  a  disjoint 


component  is  created  which  hcis  no  loop.  If  IR  and  KEY(IR)  do  not  appear  on  a  root 
loop  of  Qfu-  ^  disjoint  component  is  created,  just  as  in  the  (PNSC)  algorithm. 
If  IR  and  KEY(IR)  do  appear  on  a  root  loop  of  Qfui  *^he  resulting  retriangulated 
component  contains  a  self-loop  rather  than  a  root  loop.  The  update  of  the  subgraph 
data  structures  is  very  similar  to  those  presented  in  Chapter  7. 

Row  and  column  additions  require  the  incorporation  of  the  incoming  node  into 
the  structure  of  ^Fu  prior  to  the  column  exchange,  just  as  in  (PNSC).  The  technique 
for  (GNSC)  is  exactly  the  same  as  in  (PNSC).  Once  the  new  node  has  been  added 
to  ^  column  exchange  operation  is  performed,  with  a  logical  arc  (which  forms 
a  self-loop  on  the  new  node  being  added  to  Ofh)  being  replaced  by  the  new  arc 
(column)  being  added. 

Finally,  row  exchanges  are  handled  as  a  two-stage  update,  exactly  as  in  (PNSC). 

The  task  of  determining  the  singularity  of  Fn  ( Is_Factored_KerneljSingu- 
lar)  is  more  challenging  in  (GNSC)  than  in  either  of  the  other  two  algorithms  due 
to  the  arbitrary  nonzero  elements.  In  each  of  the  previous  two  algorithms,  we  have 
been  able  to  determine  with  certainty  the  singularity  of  Fu  indirectly.  In  the  (GUB) 
factorization,  Fn  =  An,  a  signed  identity  matrix.  Since  An  is  orthogonal  (i.e., 
Af,-AH  =  /),  F„  is  perfectly  conditioned  (see  c.g.,  Golub  and  Van  Loan  [1983]  for 
a  discussion  of  matrix  conditioning).  In  (PNSC),  Fn  may  be  triangulated,  placing  it 
in  a  form  with  plus  and  minus  ones  on  the  diagonal.  Thus,  the  determinant  of  Fn  is 
either  plus  or  minus  one.  We  may  therefore  determine  its  singularity  by  examining 
the  structure  of  Qfu  rather  than  considering  Fn  directly.  In  either  case,  as  long 
as  the  accumulated  round-off  error  in  the  current  representation  of  the  problem 
is  such  that  we  can  distinguish  plus  or  minus  one  from  zero,  we  may  rely  on  the 
structure  of  our  representation  of  to  deduce  the  singularity  of  Fn-  Thus,  we 
are  able  to  discern  singularity  logically  and  need  not  resort  to  analytic  methods. 


In  (GNSC),  however,  we  may  not  rely  exclusively  on  the  structure  of  to  reach 
conclusions  about  the  singularity  of  Fu.  It  is  not  difficult  to  invent  examples  of 
factored  kernels  corresponding  to  instances  of  ^^GNSC)  which  are  ill-conditioned. 
We  expect  ill-conditioning  of  Fu  to  lead  to  serious  numerical  difficulties,  which 
we  wish  to  avoid  if  possible.  Since  we  generally  have  freedom  in  determining  the 
structure  of  Fn  and  since  our  fundamental  rank-one  update  of  Fu  is  the  column 
exchange,  we  use  a  heuristic  based  on  column  exchanges  to  anticipate  conditioning 
problems.  Recall  in  the  column  exchange  update,  a  column  has  been  selected  for 
removal  from  Fu  and  a  candidate  column  hais  been  identified  as  its  replacement. 
We  use  the  backpath  traversal  method  described  in  Chapter  7  to  determine  if  the 
proposed  exchange  maintains  the  required  structure  of  Qfu-  not,  we  can 

conclude  that  the  exchange  renders  Fu  singular  and  we  reject  the  candidate  arc. 

As  we  traverse  the  backpath,  we  perform  calculations  which,  if  the  exchange 
does  in  fact  preserve  the  required  structure  of  will  ultimately  compute  the 
determinant  of  the  disjoint  component  in  which  the  replacement  arc  will  appear  if  the 
proposed  exchange  is  performed.  If  the  absolute  value  of  the  computed  determinant 
is  too  small  (a  decision  controlled  by  a  implementation  parameter),  we  conclude  that 
the  proposed  update  produces  a  submatrix  of  Fu  (the  submatrix  corresponding  to 
the  new  disjoint  component)  that  may  be  ill-conditioned,  and  thus  Fu  may  itself 
be  ill-conditioned.  We  therefore  reject  the  proposed  candidate  arc. 

This  approach  is  attractive  because  it  is  computationally  inexpensive  and, 
based  on  our  empirical  evidence,  it  works  well  in  practice.  It  is  of  course  merely 
a  heuristic,  since  it  is  easily  shown  (e.g.,  Golub  and  Van  Loan  [1983)]  that  the 
determinant  can  be  a  poor  indicator  of  matrix  condition.  Further,  we  do  not  actually 
compute  the  determinant  of  Fu,  but  rather  only  the  determinant  of  a  submatrix  of 
Fu¬ 
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The  remaining  update  actions,  Find_Column_to_Add(ICI)  and  Find_Col- 
umn_to_Remove(ICI)  are  treated  in  a  manner  analogous  to  that  in  (PNSC). 
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IX.  COMPUTATIONAL  RESULTS 


A.  INTRODUCTION 

The  algorithms  developed  here  have  been  implemented  and  extensively  tested 
within  the  framework  of  a  commercial-quality  optimization  system:  the  X- System 
of  Brown  and  Graves  [1975].  This  system  employs  the  Graves  mutual  primal-dual 
algorithm  in  a  variety  of  large  scale  optimization  applications,  including  linear,  non¬ 
linear,  mixed  integer  and  decomposed  models.  Although  we  report  computational 
results  only  for  linear  models,  our  factorizations  are  seamlessly  integrated  into  the 
X-Systerr  and  support  all  system  features. 

B.  TEST  PROBLEMS 

The  benchmark  test  suite  is  drawn  from  a  wide  variety  of  actual  applications. 
Table  (9.1)  provides  a  short  synopsis  of  each  model,  quoting  from  the  abstract 
where  a  reference  in  the  open  literature  is  available.  In  some  cases,  several  different 
instances  of  models  are  reported.  We  selected  these  models  because  they  provide 
a  representative  sample  of  size,  structure  and  taxonomy  of  contemporary  real-life 
applications  of  linear  programming.  For  those  models  that  are  mixed  integer,  we 
report  solution  statistics  for  the  initial  linear  programming  relaxation. 

TABLE  9.1:  Description  of  Test  Suite  Models 

•  GTE  “The  seven  Telephone  Operating  Companies  within  GTE  have  adopted 
an  integrated  business  system  called  Capital  Program  Management  System 
(CPMS)  to  guide  their  3  billion  dollar  per  year  capital  planning.  The  system 
includes  a  large  scale  mixed  integer  programming  optimization  system  that 
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optimizes  the  critical  economic  tradeoffs  between  maximizing  the  long-term 
budget  value  of  the  firm’s  equity  and  satisfying  shorter-term  financial  con¬ 
straints,  resource  limitations  and  service  objectives.  Investment  opportunities 
for  the  next  5  years  are  modeled  as  0-1  variables  with  alternative  implemen¬ 
tations  for  each.  The  objective  is  to  maximize  the  net  present  value  of  the 
capital  portfolio.  There  are  financial  constraints  on  capital,  internally  gen¬ 
erated  funds,  net  income  to  common,  and  limits  on  resources  such  as  labor 
hours,  lines  installed,  etc.  There  are  also  constraints  that  enforce  logical  re¬ 
lationships  among  opportunities  (such  as,  if  choose  A  then  must  choose  B).” 
See  Bradley  [1986]. 

•  INVEST  “Capital  allocation  and  project  selection  for  a  large  multi-national 
firm  is  modeled  as  a  two-stage  multi-year  nonlinear  capital  budgeting  problem 
with  over  40,000  integer  variables.  Innovative  modeling  yields  subproblems 
easy  to  solve,  and  optimality  is  achieved  with  a  single  iteration  of  the  nonlinear 
master  problem.”  See  Harrison,  Bradley  and  Brown  [1989].  The  instance 
reported  here  is  a  linear  program  subproblem  of  the  two-stage  model. 

•  TANKER  crude  oil  tanker  scheduling  problem  faced  by  a  major  oil  com¬ 
pany  is  solved  using  an  elastic  set  partitioning  model.  The  model  takes  into 
account  all  fleet  cost  components,  including  the  opportunity  cost  of  ship  time, 
port  and  canal  charges,  demurrage  and  bunker  fuel.  The  model  determines 
optimal  speeds  for  the  ships  and  the  best  routing  of  ballast  (empty)  legs,  as 
well  cis  which  cargoes  to  load  on  controlled  ships  and  which  to  spot  charter.  All 
feasible  schedules  are  generated,  the  cost  of  each  is  accurately  determined  and 
the  best  set  of  schedules  is  selected.”  See  Brown,  Graves  and  Ronen  [1987]. 
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•  GAS  “Natural  gas  utilities  supply  about  one  fourth  of  the  energy  needs  of 
the  United  States.  From  wellhead  to  consumer,  operations  are  governed  by  an 
astounding  diversity  of  purchase,  transport  and  storage  contract  agreements 
which  enable  complex  physical  distribution  systems  to  meet  future  demands 
no  more  predictable  than  next  year’s  weather.  Gas  is  a  highly  detailed  op¬ 
timization  model  which  utilities  use  to  plan  operations  and  to  justify  such 
plans  to  regulatory  agencies  is  developed.”  See  Avery,  Brown,  Rosenkranz 
and  Wood  [1989]. 

•  ODSM  “A  commonly  occurring  problem  in  distribution  system  design  is  the 
optimal  location  of  intermediate  distribution  facilities  betw’een  plants  and  cus¬ 
tomers.  A  multicommodity  capacitated  single-period  version  of  this  problem 
is  formulated  as  a  mixed  integer  linear  program.  .A  solution  technique  based 
on  Benders  Decomposition  is  developed.  ...  An  essentially  optimal  solution 
is  found  and  proven  with  a  surprisingly  small  number  of  Benders  cuts.”  See 
Geoffrion  and  Graves  [1974].  The  instances  reported  here  are  decomposition 
master  problems. 

•  TAM  “The  annual  decision  on  how  much  of  the  Air  Force  procurement  budget 
should  be  spent  on  the  many  different  aircraft  and  how  much  should  be  spent 
on  the  many  different  munitions  is  of  great  interest  to  many  people.  How 
the  Air  Force  staff  develons  information  to  support  the  decision  has  changed 
over  the  years.  Currently,  a  linear  program  is  being  used  by  the  Air  Force 
Center  for  Studies  and  Analysis  and  is  being  tested  by  the  Munitions  Division 
of  the  Plans  and  Operations  Directorate  (AF/XOXFM)  for  munitions  tradeoff 
analysis.  The  LP  uses  existing  data  and  estimates  on  ( 1 )  aircraft  and  munition 
effectiveness.  (2)  target  value.  (3)  attrition,  (4)  aircraft  and  munition  costs. 
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and  (5)  existing  inventories  of  aircraft  and  munitions.  Other  fjw:tors  considered 
are  weather  and  length  of  the  conflict.”  See  Might  [1987]  and  Jackson  [1989]. 

PHOENIX  ‘“The  U.  S.  Army  operates  a  fleet  of  over  7,000  hehcopters  to 
perform  combat  and  combat  support  tasks.  Although  newer,  more  technically 
advanced  helicopters  have  been  and  are  being  procured,  the  majority  of  the 
fleet  is  still  composed  of  helicopters  that  were  built  during  the  late  1960s 
for  use  in  South  Vietnam.  These  older  airframes  are  rapidly  reaching  the 
end  of  their  useful  lives  and  must  be  (i)  replaced  by  newer,  more  advanced 
designs,  (ii)  gutted  and  refit  or  replaced  by  a  combination  of  (i)  and  (ii).  Army 
force  planners  recognize  that  this  problem  can  only  be  solved  by  a  long-term 
program  and  the  commitment  of  billions  of  dollars.  Phoenix  is  a  software 
system  employing  mixed  integer  linear  programming  to  help  Army  aviation 
staff  and  officers  develop  long-  range  helicopter  modernization  plans.”  See 
Clemence,  Teufert,  Brown  and  Wood  [1988]. 

AMM04H  “A  four-commodity  transshipment  model  for  delivery  over  time  of 
military  products  from  production  and  storage  locations  to  overseas  locations 
to  support  theater  operations  is  developed.  The  model  covers  five  physical 
echelons,  including  production  plants,  storage  depots,  ports  of  embarkation, 
ports  of  debarkation  and  geographic  field  locations.  Road,  rail,  sea  and  air 
transportation  are  modeled,  and  product  demands  are  time-  phased.  Capaci- 
tation  occurs  primarily  on  sea  and  air  links,  and  as  throughput  capacities  on 
transfer  points,  requiring  replication  of  some  echelons.  The  objective  of  the 
model  is  to  minimize  deviation  from  on-time  deliveries.”  See  Staniec  [1984]. 

GK  A  weekly  multi-plant  production/inventory/transshipment  linear  pro¬ 
gram  from  a  consumer  products  industry  is  developed.  The  model  is  meant 


to  guide  weekly  processing  and  packaging  decisions.  Production  consists  of 
two  stages:  basic  products  are  produced  and  then  packaged  into  different¬ 
sized  containers  to  yield  finished  products.  Processing  lines  typically  produce 
a  subset  of  the  basic  products  and  have  limited  capacity  with  overtime  charges 
for  weekend  shifts.  Packaging  lines  for  finished  products  are  similar.  In-house 
inventory  capacity  is  limited  although  outside  storage  is  available  at  additional 
cost.  Inter-plant  shipments  are  made  by  rail  or  truck.  See  Wood  [1989]. 

C.  METHODOLOGY 

We  wish  to  evaluate  implementations  of  our  three  algorithms  on  the  basis 
of  computation  time  and  computer  memory  requirements.  Since  each  algorithm  is 
simplex-based,  the  formal  theory  of  algorithmic  complexity  provides  no  basis  for 
preferring  one  to  another,  since  in  the  worst  case  none  enjoys  a  measure  of  running 
time  that  is  polynoniial  in  the  size  of  the  problem  specification  (see  e.g.,  Garey 
and  Johnson,  [1979]  for  a  discussion  of  algorithmic  complexity  and  Klee  and  Minty 
[1972]  for  an  analysis  of  the  Simplex  method).  Thus,  we  are  led  to  consider  ‘‘typical” 
performance  by  gathering  empirical  evidence  on  the  performance  of  the  algorithms 
solving  “typical”  problems. 

We  prefer  an  implementation  that  is  both  fast  and  requires  little  computer 
memory.  Most  researchers  who  have  reported  on  implementations  of  related  al¬ 
gorithms  have  been  concerned  primarily  with  execution  speed,  and  certainly  it  is 
important.  However,  we  have  seen  in  our  algorithmic  setting  that  once  high  speed 
storage  has  been  allocated  for  the  problem  representation  and  for  the  program  code, 
all  remaining  memory  is  available  to  store  the  representation  of  the  explicit  trans¬ 
formation  kernel  and  the  factored  kernel.  If  the  solution  trajectory  is  such  that  their 
combined  size  never  exceeds  available  memory,  we  succeed  in  solving  the  problem 


If  not,  we  fail.  When  success  or  failure  depends  on  the  total  storage  requirements 
of  the  inverse  representation,  we  may  be  willing  to  sacrifice  execution  speed  in  ex¬ 
change  for  an  economical  representation  of  the  inverse.  This  is  a  classic  theme  in 
computer  science  and  software  engineering,  and  we  believe  its  importance  in  this 
context  has  been  largely  overlooked. 

We  will  be  comparing  the  performance  of  four  separate  implementations.  The 
first  is  an  unadorned  version  of  the  X-System,  which  implements  the  Graves  mutual 
primal-dual  method  as  presented  in  Chapter  2.  There  is  no  underlying  factorization. 
We  refer  to  this  implementation  as  “X”.  The  second  implementation  is  the  GUB 
factorization  presented  in  Chapter  6,  and  referred  to  as  “GUB”.  The  third  is  the 
pure  network  factorization  of  Chapter  7,  referred  to  as  “PNET”,  and  the  last  is  the 
generalized  network  factorization,  called  “GNET”,  of  Chapter  8. 

The  ideal  approach  for  this  computational  study,  would  be  to  develop  four 
equivalent  formulations  of  each  model,  each  customized  for  it.  particular  implemen¬ 
tation  with  the  goal  of  inducing  a  large  factored  row  set  of  the  appropriate  type. 
This  approach  is  a  consistent  theme  in  the  literature  dealing  with  specialized  al¬ 
gorithms  and  one  that  we  strongly  endorse.  .Alternate  formulations  of  a  model  are 
often  available,  and  it  seems  sensible  to  choose  one  that  as  strongly  as  possible 
exploits  the  strengths  of  the  solver. 

However,  all  of  the  models  used  here  are  “off-the-shelf”  in  the  sense  that  they 
were  developed  at  various  times  by  various  modelers,  and  we  cannot  afford  to  develop 
alternate  formulations.  Thus,  our  approach  is  to  preserve  a  single,  unfactored  rep¬ 
resentation  of  each  model,  and  attempt  to  identify  favorable  row  structures  through 
the  use  of  heuristics.  Our  procedure  is  ba.sed  on  the  work  of  Brown  and  Thomen 
[1980],  Brown  and  Wright  [1983]  and  Brown,  McBride  and  Wood  [1985].  The  heuris¬ 
tics  are  greedy  and  myopic  in  the  sen.'^e  that  they  initially  consider  the  entire  row 


set  of  the  problem,  and  discard  one  row  at  a  time  without  backtracking  until  the 
remaining  set  satisfies  the  desired  row  factorization.  This  may  easily  result  in  the 
confounding  or  destruction  of  structure  intended  or  perceived  by  the  modeler.  While 
our  approach  yields  interesting  and  useful  observations  about  the  implementations, 
it  is  in  some  ways  a  poor  substitute  for  the  customized  model  method. 

Table  (9.2)  tabulates  the  important  structural  information  concerning  the 
model  instances  we  will  be  solving.  The  column  headings  may  be  interpreted  as 
follows:  m  is  the  total  number  of  structural  constraints  (note  that  this  use  is  differ¬ 
ent  from  that  in  the  problem  specifications  of  Chapters  6,  7.  and  8),  n  the  number 
of  variables,  Pgub  the  number  of  GUB  rows  identified  by  the  heuristic,  ppy  the 
number  of  pure  network  rows,  pos  the  number  of  generalized  network  ro%/s  and  NZ 
the  number  of  nonzeros  in  the  technological  coefficient  matrix  of  the  model.  For 
example,  consider  the  first  model  in  Table  (9.2),  GTE.  The  structural  constraints 
contain  57,563  nonzeros,  and  the  model  consists  of  6,624  variables.  When  viewed 
a.s  an  unfactored  mutual  primal-dual  model,  it  consists  of  960  explicit  constraints 
and  0  factored  constraints.  When  viewed  as  a  GUB  factorization,  it  consists  of  909 
factored  (GUB)  constraints  and  960  —  909  =  51  explicit  constraints.  Similarly,  when 
viewed  as  a  pure  network  factorization  (PNET),  it  contains  909  factored  rows  and  51 
explicit  rows.  Finally,  when  viewed  as  a  generalized  network  factorization  (GNET). 
it  consists  of  922  factored  rows  and  38  explicit  rows. 

D.  COMPUTATIONAL  RESULTS 

We  solved  each  of  tnese  problem  instances  on  an  IBM  3033/.'\P  under  the 
VM/CMS  operating  system  using  V'S  FORTR.AN  1.4.1.  A  virtual  machine  size 
of  six  megabytes  was  used,  simply  because  it  is  the  largest  size  normally  available 
to  us.  It  is  the  nature  of  a  time-shared  system  that  measurements  of  processing 


times  are  somewhat  imprecise  due  to  system  load  factors  and  accounting  techniques 
for  system  processing  overhead.  We  have  attempted  to  mitigate  these  effects  by 
performing  our  experiments  during  periods  of  low  system  usage. 

Table  (9.3)  displays  the  solution  times  for  each  of  the  test  problem  instances  by 
each  of  the  four  implementations.  An  indicates  that  the  solver  failed  to  solve  the 
problem  instance.  Such  a  failure  occurred  because  the  storage  requirements  for  the 
explicit  transformation  kernel  representation  exceeded  the  memory  available  and  the 
solver  terminated  gracefully.  Solution  times  represent  only  the  CPU  time  required 
to  solve  the  problem,  and  exclude  the  initial  problem  input,  the  time  required  by 
the  factorization  heuristic  to  identify  the  factored  row  structure  and  the  final  output 
to  record  the  solution.  All  figures  are  in  CPU  seconds. 

The  formulations  of  three  of  the  test  problems  were  strongly  influenced  by  the 
design  of  the  target  solver:  the  TANKER  model  possesses  a  strong  GUB  structure 
since  it  contains  a  set  of  schedule  selection  constraints  for  each  ship  (i.e.,  from  a  set 
of  candidate  schedules,  select  at  most  one).  The  AMM04H  model  is  a  multicom¬ 
modity  capacitated  transshipment  problem  and  is  thus  best  suited  to  a  pure  network 
factorization.  The  PHOENTXIO  model  design  was  shaped  by  the  generalized  net¬ 
work  factorization  paradigm.  The  nature  of  the  factored  row  structures  shown  in 
Table  (9.2)  supports  this  assertion.  In  the  T.ANKER  model,  the  factored  pure  net¬ 
work  rows  are  exactly  the  same  cls  the  the  factored  GUB  rows,  and  the  heuristic 
constructs  a  generalized  network  factored  row  set  by  identifying  one  additional  row 
to  be  paired  with  each  GUB  row.  The  .AMM04H  model  may  be  viewed  as  a  GUB 
factorization  with  a  relatively  modest  GUB  set  consisting  of  the  joint  capacitation 
constraints,  or  as  a  pure  network  factorization  with  a  relatively  large  pure  network 
(PN)  factored  row  set.  In  the  PHOENTXIO  model,  the  dominant  structure  is  clearly 
the  generalized  network  row  structure. 


As  we  expect,  the  factorization  is  most  successful  when  the  model  is  wed  to  the 
solver.  Although  we  are  surprised  to  find  the  performance  of  (P.NET)  competitive 
with  (GUB)  on  the  T.\NKER  model.  (GUB)  clearly  dominates  the  X  and  GNET 
solvers.  Similarly,  (GXET)  dominates  on  the  PHOEN'IXIO  model  and  (P.XET)  on 
the  AMM04H  model.  We  would  be  disappointed  if  the  results  were  otherwise. 

We  see  evidence  in  Table  (9.3)  to  suggest  that  the  approach  of  using  heuristics 
to  automatically  identify  factored  structure  has  its  pitfalls.  In  a  number  of  problem 
instances,  although  our  heuristics  identify  significantly  larger  factored  sets  as  we 
progress  from  the  ba.ce  system  to  (G.N'ET).  we  see  little  improvement  in  computation 
times  (INVEST,  ODSMl.  FA.MS).  In  fart,  we  see  in  the  T.ANKER  model  that 
the  temptation  to  confou  ul  the  modeler's  intended  GUB  structure  by  specifying  a 
generalized  network  factorization  leads  to  disastrous  consequences,  even  though  this 
tactic  doubles  the  size  of  the  factored  row  set.  Interestingly  enough,  when  we  apply 
the  (GNET)  solver  to  the  (GUB)  factored  row  set,  we  are  able  to  solve  the  problem 
in  16  5  CPU  seconds.  This  suggests  that  the  “quality”  of  a  row  factorization  is  not 
completely  specified  by  size  alone 

We  are  encouraged  by  the  observation  that  the  transition  from  the  basic  system 
to  (GUB)  to  (PNET)  seldom  degrades  solution  times,  even  when  doing  so  yields 
little  gain  in  the  number  of  additional  factored  rows.  This  seems  to  contradict 
popular  folklore,  which  suggests  that  computation  times  worsen  as  the  transition 
from  unadorned  Simplex  to  (GUB)  to  (PNET)  is  made  unless  the  transition  is 
accompanied  by  a  substantial  incrccise  in  the  size  of  the  factored  portion  of  the 
model.  In  fact,  computational  testing  reported  by  others  is  frequently  limited  to 
models  in  which  the  number  of  explicit  rows  is  in  the  range  of  one  to  twenty  (see 
e.g..  Chen  and  Saigal  [1977],  Glover,  Karney,  Klingman  and  Russell  [1978],  Glover 
and  Klingman  [1981]).  Our  results  are  all  the  more  remarkable  given  the  the  lack  of 


guidance  from  the  modeler  for  the  “intended"  row  factorization.  W'e  note,  however, 
that  results  are  mixed  for  the  transition  from  (PNET)  to  (GNET),  and  it  seems 
clear  that  the  applicability  of  the  (GNET)  factorization  is  not  as  general  as  that  of 
(PNET). 

Finally,  we  observe  that  in  several  models  it  is  factorization  that  separates 
success  from  failure  in  solving  the  problem  with  a  given  allocation  of  computer 
resources.  This  fact  alone  may  be  reason  to  consider  this  approach  in  practice. 

Our  second  interest  with  respect  to  computation  is  in  the  memory  requirements 
of  our  algorithms.  Our  design  strategy  allocates  memory  to  the  data  structures 
which  represent  Fn  so  that  we  may  successfully  represent  the  largest  dimension  of 
factored  kernel  that  may  possibly  arise.  Limited  memory  remains  to  store  the  repre¬ 
sentation  of  -dLi',  the  explicit  transformation  kernel.  During  the  solution  process,  if 
we  encounter  a  representation  of  which  requires  more  memory  than  is  available, 
failure  occurs.  We  measure  the  size  of  the  computer  representation  of  /iL/  *^he 
amount  of  available  computer  storage  in  terms  of  the  elements  of  that  can  be 
stored.  The  number  of  bytes  per  element  varies  according  to  the  size  of  the  problem 
(this  has  to  do  with  the  FORTR.AN  data  types  L\TEGER‘2  and  1NTEGER*4).  but 
is  generally  28  bytes  per  element.  Table  (9.!)  ii'^ts  for  each  problem  instance/solver 
pair  the  maximum  size  of  the  .dR'  representation  encountered  during  the  solution, 
measured  in  units  of  number  of  elements.  .An  asterisk  (*)  indicates  that  the  number 
shown  equals  the  maximum  number  of  units  of  storage  that  were  available,  and  thus 
failure  occurred. 

We  see  that  generally  the  representation  of  the  maximum  size  of  the  explicit 
transformation  kernel  decreases  as  the  generality  of  the  factorization  increases.  Re¬ 
calling  the  definition  of  the  explicit  transformation  kernel; 


this  trend  is  as  we  would  expect.  As  the  generality  of  the  factorization  increases, 
we  expect  the  size  of  the  factored  component  to  increase  and  the  size  of  the  ex¬ 
plicit  component  to  decrease.  Each  potentially  binding  explicit  row  which  may  be 
converted  to  a  factored  ro^/  by  adopting  a  more  general  factorization  reduces  the 
dimension  of  the  resulting  representation  of  .Also,  the  density  of  the  term 

—  £n£i^'Fi2  generally  increases  as  the  complexity  of  the  structure  of  increases. 
Assuming  the  dimension  of  is  k  by  k,  the  number  of  nonzeros  in  F{[^  in  the  GUB 
factorization  is  k.  In  the  pure  network  factorization,  the  number  of  nonzeros  in  F{[^ 
may  be  as  large  as  y,  and  in  the  generalized  network  factorization,  the  number  of 
nonzeros  in  may  be  as  large  as  n^.  We  note  that  (GNET)  again  provides  several 
exceptions  to  the  general  trend  in  Table  (9.-1). 

It  is  the  dynamic  nature  of  our  factorization  algorithms  which  marks  this 
work  as  a  departure  from  previous  research.  Table  9.5  illustrates  the  significance 
of  this  point.  The  first  column  lists  the  number  of  constraints  which  are  binding 
at  optimality,  and  the  second  column  expresses  this  as  a  percentage  of  the  total 
number  of  constraints  in  the  problem  instance.  The  results  shewn  here  are  typical 
of  real-world  large-scale  models.  It  is  usually  the  case  that  many  constraints  are 
not  binding  at  optimality,  and  there  are  computational  advantages  to  be  gained  by 
exploiting  this  fact. 

Columns  3.  4  and  5  of  Table  9.5  list  for  (GUB),  (P.N'ET)  and  (GNET)  respec¬ 
tively  the  number  of  explicit  constraints  that  are  binding  at  optimality.  Since  in 
each  implementation,  binding  factored  constraints  are  handled  more  efficiently  than 
binding  explicit  constraints,  we  sec  that  the  computational  success  of  our  dynamic 
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factorization  algorithms  is  due  to  the  fact  that  even  in  large  model  instances,  we 
are  able  to  limit  our  attention  to  a  relatively  small  number  of  explicit  constraints, 
usually  on  the  order  of  a  few  hundred  or  less.  While  this  is  well  beyond  the  size  of 
previously  reported  implementations,  our  results  show  that  it  is  quite  manageable. 

It  is  useful  to  establish  a  criteria  for  comparing  and  contrasting  the  perfor¬ 
mance  of  the  algorithms  which  accounts  for  both  execution  time  and  storage  re¬ 
quirements.  Although  more  sophisticated  models  could  undoubtedly  be  developed, 
we  offer  a  simple  model  which  we  feel  captures  the  essential  features  we  wish  to 
consider.  Define 

f(sj}  =  s-t 

where  s  is  the  total  computer  storage  (measured  in  megabytes)  required  to  solve 
a  problem  instance  (including  program  code,  original  problem  data,  tableau  rep¬ 
resentation,  factored  kernel  representation  and  explicit  transformation  kernel  rep¬ 
resentation)  and  t  is  the  execution  time  (measured  in  CPU  seconds).  f{sj)  is 
monotonically  increasing  in  s  and  (,  and  in  a  crude  fashion  it  captures  the  essential 
features  of  the  way  in  which  computer  resources  are  often  marketed  commercially. 
We  will  use  f{s,  t)  as  a  measure  of  performance.  Table  (9.6)  displays  /(s,  t)  for  each 
of  the  problem  instance/solver  combinations.  indicates  that  the  solver  failed  to 
solve  the  particular  instance. 

We  observe  that  the  trend  as  the  transitions  are  made  from  base  system  to 
(PNET)  is  a  decline  in  It  is  apparent  that  (PNET)  is  a  versatile  imple¬ 

mentation.  performing  extremely  well  on  models  with  highly  favorable  structure 
(AM.M04H,  KGl,  G.ASPX.A.  GASPiN'C)  and  comparing  favorably  with  (GUB)  and 
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(GNET)  on  every  other  model.  The  performance  of  (GN'ET)  is  generally  quite  com¬ 
parable  to  that  of  (PNET).  It  appears  to  be  about  15-20%  slower  than  (PNET) 
when  it  is  used  to  solve  a  pure  network  factorization  (GASPNC,  AMM04H).  Our 
row  elimination  heuristics  are  usually  effective  in  identifying  a  favorable  factored 
structure,  but  the  computational  results  on  the  TANKER  model  clearly  illustrate 
the  limitations  of  this  approach.  (GNET)  apparently  requires  a  more  sophisticated 
and  careful  user  than  do  (GUB)  or  (PNET),  and  in  this  sense  it  is  perhaps  a  more 
specialized  algorithm.  The  evidence  shown  here  indicates  that  (PNET)  is  a  strong 
candidate  for  use  as  a  general  implementation,  and  need  not  be  viewed  as  a  highly 
specialized  implementation  suitable  only  for  rare  model  instances.  For  the  models 
studied  here,  we  have  progressed  well  beyond  the  stage  of  solving  instances  with 
only  a  handful  of  explicit  rows. 
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TABLE  9.2:  Summary  of  Problem  Suite  Dimensions 


m 

n 

GTE 

960 

6,624 

INVEST 

1,338 

11,989 

TANKER 

83 

7,598 

GAS  PN  A 

6,848 

27,884 

GAS  PN  C 

3,794 

15,362 

GAS  PN  E 

1,184 

5,102 

GAS  GN  A 

6,848 

27,884 

GAS  GN  C 

3,974 

15,362 

GK  2 

3,819 

17,841 

GK  3 

5,728 

27,493 

GK  4 

7,636 

37,139 

ODSMl 

3,023 

11,568 

ODSM2 

591 

22,211 

TAMl 

91 

389 

TAM  2 

ISO 

1,201 

TAM  3 

269 

2,883 

TAM4 

211 

1 ,327 

TAM5 

438 

10,969 

TAMS 

420 

6,104 

TAM12 

629 

17,793 

PHOENIXIO 

1,618 

6,884 

PHOEMX30 

4,305 

4,297 

AMMO  4II 

13,963 

83,497 

Pgub 

Ppn 

Pgn 

NZ 

909 

909 

917 

57,.563 

941 

1,101 

1,168 

36,829 

33 

33 

66 

30,890 

4, .345 

5,934 

5,976 

36,702 

2,6.58 

3,418 

3,420 

19,701 

4.34 

877 

883 

7,3.55 

4,484 

5,142 

5,976 

36,702 

2,664 

3,084 

3,420 

19,701 

1,265 

2,578 

2,585 

34,809 

2,295 

3,867 

3,876 

54,289 

2,428 

5,156 

5,167 

73,766 

.523 

.540 

5.58 

21, .532 

490 

490 

490 

42,827 

22 

28 

34 

1,212 

42 

54 

66 

6,869 

63 

81 

99 

21,356 

59 

77 

98 

6,954 

102 

132 

162 

93,964 

118 

154 

196 

49,376 

177 

231 

295 

164,947 

206 

220 

1,153 

13,818 

29.3 

.303 

3,605 

16,441 

1,071 

12,892 

12,892 

166,713 

TABLE  9.3:  Solution  Times  in  CPU  Seconds 


X  GUB  PNET  GNET 


GTE 

48.4 

33.0 

37.2 

38.3 

INVEST 

2E6 

17.1 

17.1 

17.0 

TANKER 

141.8 

20.5 

17.7 

43.8 

GAS  PN  A 

* 

* 

413.2 

527.7 

GAS  PN  C 

56.9 

53.7 

35.8 

31.3 

GAS  PN  E 

4c 

146.8 

17.4 

20.7 

GAS  GN  A 

* 

♦ 

* 

521.6 

GAS  GN  C 

52.9 

50.9 

37.6 

33.7 

GK  2 

27.1 

22.8 

21.3 

19.4 

GK  3 

62.4 

58.4 

•50.4 

45.1 

GK  4 

* 

380.9 

182.2 

186.6 

Ol»m\11 

8.7 

12.1 

9.3 

8.8 

ODSM2 

34.1 

.32.5 

32.8 

30.4 

TA.Ml 

0.7 

0.7 

0.7 

0.7 

TA.M2 

2.5 

1.7 

2.2 

2.1 

TAM  3 

17.8 

14.1 

14.3 

13.5 

TAMl 

1.2 

1.2 

1.5 

1.7 

TAM5 

317.1 

300.9 

307.8 

269.6 

TAMS 

93.3 

81.8 

83.8 

77.4 

TAM  12 

6G0.9 

660.6 

624.1 

480.9 

PIIOENIXIO 

6S.6 

43.7 

30.4 

10.7 

PHOENIX30 

* 

* 

* 

495.1 

nMMO  411 

NA 

t  1,241.3 

253.9 

288.4 

•  NA  indicates  problem  instance  not  run. 

t  This  problem  was  run  on  an  IBM  -3081 K  under  the  M\  S  operating  system  using 
VS  FORTRAN  l.t.l  in  a  32  megabyte  virtual  macliine  using  an  advanced  starting 
solution.  I'lie  solution  time  shown  is  adjusted  to  account  for  the  approximate  dif¬ 
ference  in  com[)ut!!)g  speed  between  the  IBM  3()81K  and  the  IBM  .}0.).]/AI. 


TABLE  9.4:  Number  of  Elements  in  Explicit  Transformation  Kernel  Rep¬ 
resentation  of  Optimality  or  at  Failure 


X 

GUB 

PNEl 

GNET 

GTE 

11,571 

221 

214 

196 

INVEST 

8,380 

682 

654 

521 

TANKER 

2,256 

415 

312 

81 

GAS  PN  A 

*132,995 

*130,883 

35,395 

65,558 

GAS  PN  C 

21,912 

14,765 

260 

364 

GAS  PN  E 

*174,482 

111,821 

4.808 

4,312 

GAS  GN  A 

*132,995 

*130,883 

*127,780 

65,081 

GAS  GN  C 

17,966 

8.017 

1,358 

608 

GK  2 

8,280 

2,719 

221 

184 

GK  3 

13,991 

7.771 

388 

310 

GK  4 

*110,334 

79,646 

6,532 

9,155 

ODSMl 

1,322 

861 

259 

421 

ODSM2 

418 

3 

3 

3 

TAMl 

1.713 

1.231 

841 

1,165 

TAM  2 

3,042 

2.078 

1.444 

1.610 

TAM3 

10,430 

8,326 

5,305 

6,525 

TAM4 

1,489 

1,302 

867 

1 .332 

TAMS 

29.661 

23,790 

15.689 

18,737 

TAM8 

19,545 

14,886 

8,351 

11,182 

TAM12 

37.716 

28,935 

16,097 

20,486 

PHOENIXIO 

73,737 

38,705 

21.879 

1,4  58 

PHOENIX30 

*153,781 

*151,746 

*148,348 

11,523 

AMMO  4H 

NA 

122,085 

23 

23 

TABLE  9.5:  Summary  of  the  Number  of  Binding  Explicit  Constraints  at 
Optimality 


Total 

Binding 

Percent 

GUB 

PNET 

GNET 

GTE 

552 

57.5 

19 

18 

18 

INVEST 

758 

.56.7 

199 

194 

162 

TANKER 

50 

60.2 

31 

30 

9 

GAS  PN  A 

3,051 

44.6 

292 

319 

GAS  PN  C 

2,337 

61.6 

849 

66 

91 

GAS  PN  E 

897 

75.8 

572 

88 

81 

GAS  GN  A 

3,045 

44.5 

4c 

3.52 

GAS  GN  C 

2,331 

.58.7 

863 

402 

88 

GK  2 

1 ,950 

.51.1 

1,011 

93 

90 

GK  3 

2,912 

.50.8 

1 ,360 

140 

141 

GK  1 

4,119 

53.9 

2,477 

363 

3.56 

ODSMl 

297 

9.8 

53 

49 

47 

0DS.M2 

448 

75.4 

219 

215 

0 

TAMl 

46 

.50.5 

44 

32 

35 

TA.\12 

87 

48.3 

77 

51 

61 

TAM3 

148 

.55.0 

131 

97 

107 

TAM4 

121 

.57.3 

110 

59 

81 

TAM  5 

252 

57.5 

228 

170 

186 

TAMS 

276 

65.7 

238 

140 

174 

TAM  12 

413 

65.7 

355 

208 

252 

PHOENIXK) 

1 ,085 

67.1 

1,083 

1,082 

77 

P11OEN1X30 

3,477 

80.8 

♦ 

* 

109 

AMMO  4H 

2,889 

20.7 

2,882 

7 

7 

l(.l 


TABLE  9.6:  f(s,t)  in  Megabyte  seconds 


X 


GTE 

84 

INVEST 

25 

TANKER 

II2 

GAS  PN  A 

4c 

GAS  PN  C 

73 

GAS  PN  E 

* 

GAS  GN  A 

* 

GAS  GN  C 

62 

GK  2 

27 

GK  3 

90 

GK  4 

* 

ODSMl 

4 

OI)SM2 

49 

TAMl 

0.1 

TAM2 

1 

TAM3 

14 

TAM  4 

0.1 

TAM5 

750 

TAMS 

133 

TAM12 

2,376 

PHOENIXIO 

169 

PHOENTX30 

* 

AMMO  4H 

NA 

GUB 

PNET 

GNET 

47 

54 

64 

16 

19 

20 

15 

13 

42 

* 

926 

1 ,77 1 

58 

27 

31 

509 

9 

14 

4c 

* 

1,774 

45 

30 

34 

19 

18 

21 

74 

60 

66 

1,368 

309 

387 

6 

5 

7 

47 

48 

53 

0.1 

0.1 

0.1 

1 

1 

1 

10 

9 

12 

0.1 

0.1 

0.1 

662 

614 

622 

106 

95 

110 

2,213 

1,882 

1,632 

65 

32 

7 

4c 

4c 

693 

9,093 

933 

1194 

X.  CONCLUSIONS 


We  have  presented  three  dynamic  row  factorization  algorithms  for  solving 
large-scale  linear  programs.  Although  each  may  be  used  to  solve  any  LP  instance, 
each  is  designed  to  exploit  a  particular  model  row  structure:  generalized  upper 
bounds  (GUB),  pure  network  rows  (PNET)  or  generalized  network  rows  (GNET). 

Previous  research  by  others  generally  suggests  that  specialized  algorithms  such 
as  those  presented  here  are  useful  only  when  the  factored  structure  completely  dom¬ 
inates  the  structure  of  the  model  instance.  There  are  reports  of  algorithms  for 
solving  problems  having  a  single  unfactored  (explicit)  constraint  (Hultz  and  Kling- 
man  [1978],  Klingman  and  Russell  [1978]).  When  implementations  are  reported, 
problem  suites  are  limited  to  instances  having  a  very  small  number  of  explicit  con¬ 
straints,  typically  in  the  range  from  one  to  twenty  (Chen  and  Saiga!  [1977],  Glover, 
Karney,  Klingman  and  Russell  [1978],  Glover  and  Klingman  [1981]).  The  consensus 
seems  to  be  that  such  algorithms  are  appropriately  viewed  as  specialized  algorithms, 
useful  only  for  solving  very  special  problem  instances. 

Our  experience  strongly  refutes  this  view.  We  find  the  performances  of  our  im¬ 
plementations  of  the  dynamic  factorization  algorithms  are  competitive  with  that  of 
a  commercial-quality  optimization  system  on  every  model  instance  we  have  tested. 
This  is  particularly  remarkable  for  two  rea.sons.  First,  our  test  suite  consists  of 
models  developed  by  skilled  modelers  specifically  to  exploit  the  capabilities  and 
characteristics  of  the  solver  with  which  our  implementations  are  competing.  Sec¬ 
ond.  we  must  select  the  row  factorizations  without  the  benefit  of  guidance  from  the 
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modeler,  relying  instead  on  useful  but  imperfect  heuristics.  Despite  these  computa¬ 
tional  handicaps,  our  tests  show  our  implementations  to  be  at  least  as  efficient  as  a 
well-respected  commercial-quality  optimization  system. 

Our  development  has  stressed  the  similarity  between  the  algorithms  and  the 
natural  extension  which  leads  from  one  to  the  next.  This  is  in  contrast  to  the  devel¬ 
opment  that  has  been  reported  for  similar,  non-dynamic  algorithms  (e.g.,  Dantzig 
and  Van  Slyke  [1967],  Klingman  and  Russell  [1978]  and  Hultz  and  Klingman  [1978]) 
in  which  the  specifics  of  the  individual  algorithm  obscure  the  generality  of  the  ap¬ 
proach.  The  conceptual  difference  between  our  algorithms  is  seen  to  be  largely 
isolated  to  the  structure  of  a  single  algebraic  entity,  the  factored  kernel.  By  abstract¬ 
ing  the  structure  of  the  factored  kernel  and  concentrating  on  the  general  algorithm 
design,  we  demonstrate  the  versatility  and  flexibility  of  this  approach. 

We  are  gratified  to  find  that  the  modularity  suggested  by  the  algorithmic 
development  can  be  realized  in  an  implementation  design.  We  succeed  in  developing 
a  software  suite  which  displays  a  “single-system  image”.  The  modularity  of  the 
algorithm  allows  the  definition  of  an  “abstract  data  type”  (see,  e.g.,  Aho,  Hopcroft 
and  Ullman  [1974])  which  isolates  the  data  structures  and  update  procedures  for  the 
factored  kernel  from  the  rest  of  the  implementation.  Each  factorization  is  seamlessly 
integrated  within  the  system  design,  presenting  a  single  design  image. 

The  early  1980's  produced  a  great  deal  of  r^^search  in  the  area  of  automatic 
identification  of  special  structure  in  LP  models  (see,  e.g.,  Gunawardane  and  Schragc 
[1977],  Glover  [1980].  Schrage  [1981],  Brown.  McBride  and  Wood  [1985]  and  Bixby 
and  Fourer  [1986]).  We  ha^■e  incorporated  the  most  useful  of  these  ideas  into  our 
implementation,  and  we  have  what  we  believe  to  be  the  first  implementation  which 
supports  the  automatic  identification  of  factored  row  sets.  This  capability  may 
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be  used  to  identify  new  factored  structure  or  to  validate  or  augment  a  modeler- 
provided  recommendation.  Our  computational  experience  indicates  that  while  this 
approach  is  not  a.s  promising  as  perhaps  first  envisioned,  it  is  nonetheless  a  valuable 
tool.  When  faced  with  the  choice  of  either  solving  an  unfactored  model  instance 
or  automatically  identifying  a  factored  structure  and  then  using  the  corresponding 
solver,  our  results  show  that  the  latter  is  nearly  always  to  be  preferred.  Our  results 
seem  to  suggest,  however,  that  in  addition  to  quantity  of  factored  rows,  the  issue 
of  quality  of  factored  rows  exerts  influence  on  the  performance  of  the  factorization 
algorithms.  While  not  well  understood,  it  is  clear  that  the  myopic  approach  of 
our  heuristics  is  no  substitute  for  the  modeler’s  guidance  in  identifying  factored 
structure. 

Several  areas  suggest  themselves  for  further  research.  Certainly  additional  fac¬ 
tored  structures  can  be  examined.  For  example,  one  approach  for  treating  factored 
column  structures  (“complicating  columns")  is  to  allow  the  partitioning  of  rows  into 
the  categories  “factored”  and  “explicit"  to  vary  as  the  algorithm  progresses.  That 
is,  allow  factored  row  set  membership  to  be  determined  with  respect  to  the  column 
structure  of  the  currently  nonbasic  variables  rather  than  \’itf  respect  to  the  column 
structure  of  all  problem  variables.  While  conceptually  simple,  such  a  generalization 
seems  to  present  significant  algorithmic  challenges. 

General  algorithms  are  sometimes  useful  in  specialized  contexts.  For  example, 
processing  networks  (Koene  [1982])  are  network  models  which  allow  proportional 
flow  restrictions  on  the  arcs  entering  or  leaving  some  nodes.  One  formulation  of  such 
a  model  results  in  a  pure  or  generalized  network  structure  with  a  set  of  complicating 
columns.  Chen  and  Engquist  [1986]  propose  a  primal  partitioning  algorithm  for 
solving  processing  network  problems.  An  alternate  formulation  yields  a  pure  or 
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generalized  structure  with  complicating  rows,  and  we  note  that  this  iS  precisely  the 
structure  we  seek  for  our  network  factorizations. 

The  multicommodity  capacitated  transshipment  problem  (MCTP)  has  been 
the  subject  of  much  research  over  the  years,  and  a  number  of  specialized  algorithms 
(see,  e.g.,  Assad  [1978]  or  Kennington  [1978])  have  been  proposed  to  solve  it.  Adopt¬ 
ing  a  general  perspective,  MCTP  may  be  viewed  either  as  a  GUB  model  or  as  a  pure 
network  model  with  side  constraints,  and  either  view  might  be  preferred  depending 
upon  the  dimensions  of  the  particular  instance  under  consideration.  Our  compu¬ 
tational  experience  indicates  that  the  pure  network  factorization  algorithm  offers  a 
powerful  technique  for  solving  MCTP.  As  an  experiment,  we  customized  our  (PN) 
implementation  to  exploit  the  special  structure  of  the  side  constraints  in  MCTP.  It 
is  interesting  to  note  that  in  our  scientific  computing  environment,  we  observed  no 
difference  in  solution  times  between  the  customized  version  of  (PN)  and  the  original 
implementation. 

Finally,  all  the  approaches  we  have  considered  assume  the  prior  existence  of  a 
specific  structure  in  the  factored  rows  which  in  turn  determines  the  structure  of  the 
factored  kernel.  An  extension  of  this  general  approach  is  to  relax  the  requirement  for 
strict  conformance  to  a  specific  structure.  Instead,  we  might  allow  the  factored  row 
structure  to  be  “nearly”  homogeneous.  For  example,  we  may  allow  a  small  number 
of  complicating  columns  to  disturb  what  is  otherwise  a  factored  pure  network  row 
structure.  VV'e  then  expect  the  structure  of  the  factored  kernel  to  be  dominated  by 
that  induced  by  the  predominant  row  structure,  with  only  occasional  complications 
due  to  the  excep.ional  row  structure.  We  allow  for  this  exceptional  structure  in 
the  factored  kernel  by  identifying  it  “on-the-fly'  as  the  algorithm  progresses,  and 
treating  it  in  an  appropriate  manner.  This  approach  may  be  thought  of  as  a  hy¬ 
brid  between  the  factored  method  developed  here  and  dynamic  ba'^is  triangulation 


methods  (see,  e.g..  Hellermap  and  Rarick  [1971]  and  [1972],  Saunders  [1976]  and 
McBride  [1980]). 

Dynamic  extrinsic  factorization  is  subsumed  if  we  activate  functions  in  the 
update  analogous  to  the  secondary  exchanges  now  employed.  Essentially  all  that 
has  to  be  done  is  ensure  that  successive  factored  components  retain  their  stipulated 
special  structure.  In  our  estimation,  this  will  only  be  justified  in  cases  where  the 
model  structure  is  amenable,  and  quite  likely  will  require  some  model-specific  fea¬ 
tures  to  perform  well  on  difficult  models.  We  have  limited  our  experimentation  to 
those  static  extrinsic  cases  which  are  believed  to  be  most  useful. 
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